Graph-Guided Banding of the Covariance Matrix

arXiv:1606.00451v2 [stat.ME] 15 Feb 2018

Jacob Bien∗
Data Sciences and Operations
Marshall School of Business
University of Southern California

Abstract
Regularization has become a primary tool for developing reliable estimators of the covariance
matrix in high-dimensional settings. To curb the curse of dimensionality, numerous methods
assume that the population covariance (or inverse covariance) matrix is sparse, while making
no particular structural assumptions on the desired pattern of sparsity. A highly-related, yet
complementary, literature studies the specific setting in which the measured variables have a
known ordering, in which case a banded population matrix is often assumed. While the banded
approach is conceptually and computationally easier than asking for “patternless sparsity,” it
is only applicable in very specific situations (such as when data are measured over time or
one-dimensional space). This work proposes a generalization of the notion of bandedness that
greatly expands the range of problems in which banded estimators apply.
We develop convex regularizers occupying the broad middle ground between the former
approach of “patternless sparsity” and the latter reliance on having a known ordering. Our
framework defines bandedness with respect to a known graph on the measured variables. Such
a graph is available in diverse situations, and we provide a theoretical, computational, and
applied treatment of two new estimators. An R package, called ggb, implements these new
methods.

Keywords: covariance estimation, high-dimensional, hierarchical group lasso, network

1

Introduction

Understanding the relationships among large numbers of variables is a goal shared across many scientific areas. Estimating the covariance matrix is perhaps the most basic step toward understanding
these relationships, and yet even this task is far from simple when sample sizes are small relative to
the number of parameters to be estimated. High-dimensional covariance estimation is therefore an
active area of research. Papers in this literature generally describe a set of structural assumptions
that effectively reduces the dimension of the parameter space to make reliable estimation tractable
given the available sample size. Some papers focus on eigenvalue-related structures and shrinkage
(e.g., Johnstone 2001, Fan et al. 2008, Won et al. 2013, Donoho et al. 2013, Ledoit & Wolf 2004,
among many others); others assume that correlations are constant or blockwise constant (such as
in random effects models); and a very large number of papers introduce sparsity assumptions on
∗
The author thanks Christian Müller for a useful conversation and acknowledges the support of the NSF grant
DMS-1405746 and NSF CAREER Award DMS-1748166.

1

the covariance matrix (or its inverse). Still other papers suggest combinations of these structures
(Chandrasekaran et al. 2010, Bien & Tibshirani 2011, Luo 2011, Rothman 2012, Fan et al. 2013,
Liu et al. 2013, Xue et al. 2012).
Given the large number of papers focused on the sparsity assumption, it is perhaps surprising
to note that nearly all such papers fall into one of just two categories:
1. Methods that place no assumption on the pattern of sparsity but merely require that
the number of nonzeros be small. This assumption is placed on either the inverse covariance
matrix (Dempster 1972, Meinshausen & Bühlmann 2006, Yuan & Lin 2007, Rothman et al.
2008, Banerjee et al. 2008, Friedman et al. 2008, Peng et al. 2009, Yuan 2010, Cai et al. 2011,
Khare et al. 2014) or on the covariance matrix itself (Bickel & Levina 2008b, Rothman et al.
2009, Lam & Fan 2009, Bien & Tibshirani 2011, Rothman 2012, Xue et al. 2012, Liu et al.
2013).
2. Methods that assume the variables have a known ordering. These methods generally
assume a banded (or tapered) structure for the covariance matrix, its inverse, or the Cholesky
factor (Wu & Pourahmadi 2003, Huang et al. 2006, Levina et al. 2008, Bickel & Levina 2008a,
Wu & Pourahmadi 2009, Rothman et al. 2010, Cai et al. 2010, Cai & Yuan 2012, Bien et al.
2015a, Yu & Bien 2017).
We can view these two categories of sparsity as extremes on a spectrum. Methods in the first
category tackle a highly challenging problem that is widely applicable. Being completely neutral
topthe
pattern of zeros leads to an enormous number of sparsity patterns to choose from: there are
(2) for an s-sparse model. By constrast, methods in the second category address a more modest
s
problem that is much less widely applicable. Using a simple banding estimator, there are only p
sparsity patterns to choose from (one for each bandwidth) and there are few examples other than
data collected over time in which such estimators apply.
This paper introduces a new category of sparsity patterns that aims to bridge the wide divide between these two categories of sparsity pattern. The goal is to produce an estimator that
combines the problem-specific, targeted perspective of banded covariance estimation with the wide
applicability of the patternless sparsity approaches. The key idea is to generalize the notion of
banding beyond situations in which the variables have a known ordering. In particular, we replace
the assumption that the variables “have a known ordering” with the assumption that they “lie on
a known graph.” The latter is a generalization because variables having a known ordering can be
expressed as variables lying on a path graph.
To be precise, we will assume that we observe n independent copies of a random vector X ∈ Rp ,
where the p variables are associated with a known graph G = ([p], E). We will refer to G as the seed
graph. The right panel of Figure 1 shows a traditional B-banded matrix. The term “banded” refers
to the diagonal strip of nonzero elements that are within some distance B (called the “bandwidth”)
of the main diagonal. The middle panel of Figure 1 shows the graph with this adjacency matrix.
Observing that this graph, denoted GB , is the B-th power of the path graph shown in the left
panel of Figure 1 suggests the following generalization of the term “banded.” We will use dG (j, k)
to denote the length of the shortest path between nodes j and k on G.
Definition 1. We say that a matrix Σ is B-banded with respect to a graph G if supp(Σ) ⊆
E(GB ), that is Σjk = 0 if dG (j, k) > B.

2

●

●

●
●

●

●

●

●

●

Seed graph

●

●

●

●
●

●

●

●

●

●

●

●

5

●

B=4

●

●

●

10

●
15

●

●
●

●

●

●
●

●
20

●

●

●

●

●

●

5

10

15

20

Figure 1: A matrix with bandwidth B = 4 (right panel) can be thought of as the adjacency matrix of the
B-th power (center panel) of a path graph G (left panel). To form the B-th power of a graph G, one connects
any pair of nodes that are within B hops of each other on G.
●

●
●

●

●

●

●
●

●
●

●
●

5

●

Seed graph

B=4

●

10

●
●

●

●
●

●

●

●

●

●

●
●

●

●

●
●

●
●

15

●

●

●
●

●

20

●

●

●

5

10

15

20

Figure 2: Given a seed graph G (left panel), one forms GB (B = 4 case shown in center panel) by connecting
any pair of nodes within B hops of each other in G. Definition 1 refers to the resulting adjacency matrix
(right panel) as B-banded with respect to the graph G.

Bickel & Lindner (2012) have a similar notion, for a metric ρ, which they call ρ-generalized
banded operators.
Figure 2 shows the sparsity pattern of a 4-banded matrix with respect to a seed graph. The
bandwidth B can vary from B = 0, which corresponds to Σ being a diagonal matrix, to B being
the diameter of G, which corresponds to Σ being a dense matrix when G is a connected graph
(when G is not connected, the densest Σ attainable is, up to permutation of the variables, a blockdiagonal matrix with completely dense blocks). Despite the simple interpretation of taking powers
of a seed graph, the resulting sparsity patterns can be quite complex and varied depending on the
seed graph. There are many applications in which measured variables lie on a known graph that
would be relevant to incorporate when estimating the covariance matrix Σ:
• Image data: Consider a p1 × p2 pixeled image. The p = p1 p2 pixels lie on a two-dimensional
lattice graph. Nearby pixels would be expected to be correlated.
• fMRI data: One observes three-dimensional images of the brain over time. The voxels (threedimensional pixels) lie on a three-dimensional lattice.
3

●

●

●

●
●

●
●

●
●

B5 = 1
B8 = 1
B6 = 4
B4 = 2
B3 = 2

●
●

●

●

Seed graph

●

5

10

●
●

●

●
●

●

●

●

●

●

●
●

●

●

●
●

●

15

●
●

●

●
●

●

20

●

●

●

5

10

15

20

Figure 3: Given a seed graph G (left panel) and a set of node-specific bandwidths, one forms a
new graph (shown in center panel) by connecting any pair of nodes within Bj hops of node j in
G (five of the 20 nodes have nonzero bandwidths with values shown in figure; the remaining 15
nodes have zero bandwidths). Definition 2 refers to the resulting adjacency matrix (right panel) as
(B1 , . . . , Bp )-banded with respect to the graph G.
• Social network data: A social media company wishes to estimate correlations between users’
purchasing behaviors (e.g., do they tend to buy the same items?). The social network of users
is the known seed graph.
• Gene expression data: Genes that encode proteins that are known to interact may potentially
be more likely to be co-expressed than genes whose proteins do not interact. Biologists have
constructed protein interaction networks (Das & Yu 2012), which may be taken as seed
graphs, leveraging protein interaction data for better estimation of the covariance matrix of
gene expression data.
In each example described above, the naive approach would be to vectorize the data and ignore
the graph structure when estimating the covariance matrix. However, doing so when the truth is
banded with respect to the known graph would be unnecessarily working in the patternless sparsity
category when in fact a much simpler approach would work.
In some situations, a known seed graph is available, but the assumption of a single bandwidth B
applying globally is unrealistic. We can therefore extend the definition of bandedness with respect
to a graph as follows:
Definition 2. We say that a matrix Σ is (B1 , . . . , Bp )-banded with respect to a graph G, if
Σjk = 0 for all (j, k) satisfying dG (j, k) > max{Bj , Bk }.
In words, this definition says that each variable has its own bandwidth, which defines a neighborhood size, and it cannot be correlated with any variables beyond this neighborhood. For example,
on a social network, the bandwidth can be thought informally as having to do with the influence of
a person. A very influential friend might be correlated with not just his friends but with his friends’
friends (this would correspond to Bj = 2). By contrast, a “lone wolf” might be uncorrelated even
with his own friends (Bj = 0). Figure 3 shows an example of a sparsity pattern emerging from the
same seed graph as in Figure 2 but with variable bandwidths.
We propose two estimators in this paper that make use of a latent overlapping group lasso
penalty (Jacob et al. 2009, Obozinski et al. 2011), which is used in situations in which one desires
4

a sparsity pattern that is a union of a set of predefined groups. Both Definitions 1 and 2 can
be described as such, so the difference between the two methods lies in the choice of groups. In
Section 2, we introduce and study the graph-guided banding estimator, which produces estimates of
the covariance matrix that are banded in the sense of Definition 1. The limited number of sparsity
patterns attainable by this estimator makes it potentially undersirable in certain situations, and
we therefore introduce and study (in Section 3) the graph-guided banding with local bandwidths
estimator, which attains bandedness in the more flexible sense of Definition 2. Section 4 compares
these methods empirically both in simulation and in naturally occurring data. Section 2.3 concludes
with a discussion of related work, insights learned, and future directions.
To our knowledge, there is little work on using known graphs in high-dimensional covariance
estimation. Cai & Yuan (2015) study minimax lower bounds and adaptive estimation of covariance matrices when variables lie on a d-dimensional grid. Deng & Yuan (2009) consider inverse
covariance estimation in the setting where variables lie in a metric space. Their approach uses a
non-negative Garrotte semidefinite program formulation to perform inverse covariance estimation
in which inverse covariance decays with increasing distance between variables. The method requires
an initial estimate of the inverse covariance matrix (when n > p, one takes S −1 ) and then shrinks
so that variables farther apart have smaller partial covariances.
Notation: Given a p × p matrix A and a set of indices g ∈ {1, . . . , p}2 , we write Ag to denote the
p × p matrix that agrees with A on g and is zero on g c , and we write |g| for the cardinality of the
set g. The graph GB is the Bth power of G, meaning the graph in which any nodes that are within
B of each other in G are connected. The constraint Σ  δIp means that the smallest eigenvalue
of Σ is at least δ. Given a square matrix Σ, the matrix Σ− = Σ − diag(Σ) agrees with Σ on
off-diagonals but is zero on the diagonal. We write an . bn to mean that there exists a constant
C > 0 for which an ≤ Cbn for n sufficiently large.

2

Graph-guided banding with global bandwidth

Definition 1 introduces the simplest notion of bandedness with respect to a graph. In particular,
a single bandwidth B is postulated to apply globally across the entire graph. This means that the
support of Σ is given by one of the following groups of parameters:
gb = {jk : 1 ≤ dG (j, k) ≤ b}

for b = 1, . . . , diam(G),

(1)

where diam(G) is the diameter of the seed graph G. The group lasso (Yuan & Lin 2006) is a
well-known extension of the lasso (Tibshirani 1996) that is useful for designing estimators in which
groups of variables are simultaneously set to zero. Given a sample
of n independent random vectors
P
X1 , . . . , Xn ∈ Rp , form the sample covariance matrix, S = ni=1 (Xi − X̄n )(Xi − X̄n )T . We define
the graph-guided banding estimator with global bandwidth (GGB-global) as


1
2
− c
Σ̂λ = arg min
kS − ΣkF + λP (Σ; G) s.t. [Σ ]gM = 0 ,
(2)
ΣδIp 2
where
P (Σ; G) =

min
{V(b) ∈Rp×p }

(M
X

wb kV(b) kF

b=1

s.t. ΣgM =

M
X
b=1

5

)
Vg(b)
b

is the latent overlapping group lasso penalty (Jacob et al. 2009) with group structure given by (1).
(b)
Each Vgb matrix has the “graph-guided” sparsity pattern gb , so heuristically, this penalty favors
Σ̂λ that can be written as the sum of a small number of such sparse matrices. While the penalty
is itself defined by an optimization problem, in practice we can solve (2) directly without needing
to evaluate P (Σ; G), as described in the next section.
The parameter M will typically be taken to be diam(G); however, in certain applications in
which diam(G) is very large and a reasonable upper bound on B is available, it may be compuc = 0 only has an effect
tationally expedient to choose a much smaller M . The constraint [Σ− ]gM
when M < diam(G). Of course, one must be cautious not to inadvertently choose M < B since
then the resulting estimates will always be sparser than the true covariance matrix. Thus, if one
does not have a reliable upper bound for B, one should simply take M = diam(G). The tuning
parameter λ ≥ 0 controls the sparsity level of Σ̂λ (with large values of λ corresponding to greater
sparsity). The wg ≥ 0 are weights that control the relative strength of the individual `2 norms
in the penalty. In studying the theoretical properties of the estimator (Section 2.2), it becomes
clear that wb = |gb |1/2 is a good choice. The parameter δ is a user-specified lower bound on the
eigenvalues of Σ̂λ . While our algorithm is valid for any value of δ, the most common choice is
δ = 0, which ensures that Σ̂λ is a positive semidefinite matrix. In certain applications, a positive
definite matrix may be required, in which case the user may specify a small positive value for δ.
In the special case that G is a path graph (and dropping the eigenvalue contraint), the GGB-global estimator reduces to the estimator introduced in Yan & Bien (2017). In the next section, we show
that (2) can be efficiently computed.

2.1

Computation

The challenge of solving (2) lies in the combination of a nondifferentiable penalty and a positive
semidefinite constraint. Without the penalty term, Σ̂λ could be computed with a single thresholding
of the eigenvalues; without the eigenvalue constraint, Σ̂λ is nothing more than the proximal operator
of λP (·; G) evaluated at S. The following lemma shows that we can solve (2) by alternate application
of the proximal operator of λP (·; G) and a thresholding of the eigenvalues. We present this result
for general convex penalties since the same result will be used for the estimator introduced in
Section 3.
Lemma 1. Consider the problem

Σ̂ = arg min

ΣδIp


1
2
kS − ΣkF + λΩ(Σ) ,
2

(3)

where Ω is a closed, proper convex function, and let


1
2
ProxλΩ (M) = arg min
kM − ΣkF + λΩ(Σ)
Σ
2
denote the proximal operator of λΩ(·) evaluated at a matrix M. Algorithm 1 converges to Σ̂.
Proof. See Appendix A.
In Appendix A, we show that Algorithm 1 corresponds to blockwise coordinate ascent on the
dual of (3). While the alternating direction method of multipliers (ADMM) would be an alternative
6

Algorithm 1 Algorithm for solving (3)
• Initialize C.
• Repeat until convergence:
1. B ← (S + C) − ProxλΩ (S + C)
2. C ← U · diag([Λii + δ]+ ) · UT where B − S = UΛUT is the eigendecomposition.
• Return Σ̂ ← S − B + C.

approach to this problem, we prefer Algorithm 1 because, unlike ADMM, its use does not require
choosing a suitable algorithm-specific parameter. Invoking the lemma, all that remains for solving
(2) is to specify an algorithm for evaluating ProxλP (·;G) . The group structure (1) used in this
estimator is special in that
g1 ⊂ g2 ⊂ · · · ⊂ gM .
Yan & Bien (2017) study the use of the latent overlapping group lasso with such hierarchically
nested group structures and show that, for group structures of this form, the proximal operator
can be evaluated essentially in closed form (which is not true more generally). Algorithm 3 of Yan
& Bien (2017) shows that evaluating the proximal operator amounts to identifying a sequence of
“breakpoints” k0 < k1 < . . . < km and then applying groupwise soft-thresholding to the elements
in gki+1 − gki for i = 0, . . . , m − 1.
In practice, we solve (3) for a decreasing sequence of values of λ. We initialize C in Algorithm
1 as 0 for the first value of λ and then for each subsequent value of λ we use the value of C from
the previous λ as a warm start.

2.2

Theory

In this section, we study the statistical properties of the GGB estimator (2). For simplicity, we will
focus on the special case in which δ = −∞ (or, equivalently, we drop the eigenvalue constraint);
however, in Section 2.2.3 we provide conditions under which our results apply to the δ = 0 case.
We introduce two conditions under which we prove these results.
Assumption 1 (distribution of data). We draw n i.i.d. copies of X ∈ Rp , having mean 0 and
covariance Σ∗ . The random vector has bounded variances, maxj |Σ∗jj | ≤ κ̄, and sub-Gaussian
√ ∗
2
marginals: E[etXj / Σjj ] ≤ eCt for all t ≥ 0. Both C and κ̄ are constants.
Assumption 2 (relation between n and p). We require that p ≤ eγ1 n for some γ1 > 0.
Assumptions 1 and 2 are sufficient for establishing that the sample covariance matrix S concentrates
around Σ∗ with high probability for large n. These assumptions are made in Bien et al. (2015a),
and similar assumptions appear throughout the literature.
2.2.1

Frobenius norm error

We begin by producing the rate of convergence of our estimator in Frobenius norm.

7

Theorem 1. Suppose Σ∗ is B ∗ -banded with q
respect to the graph G (in the sense of Definition 1)
and Assumptions 1–2 hold. Then with λ = x log(max{p,n})
(for x sufficiently large), wb = |gb |1/2 ,
n
δ = −∞, and M = diam(G), there exists a constant c > 0 such that
kΣ̂λ − Σ∗ k2F .

(p + |gB ∗ |) log(max{p, n})
n

with probability at least 1 − c/ max{p, n}.
Proof. See Appendix B.
The term p + |gB ∗ | is the sparsity level of Σ∗ (with the initial p representing the diagonal
elements). Theorem 1 shows that when p < n and moreover p + |gB ∗ | is o(n/ log n), the GGB
estimator is consistent in Frobenius norm. By contrast, consistency of the sample covariance
matrix S in this setting requires that p2 be o(n/ log n). Thus, when |gB ∗ |  p2 , GGB is expected
to substantially outperform the sample covariance matrix. When p > n, the theorem does not
establish consistency of GGB. This fact is not, however, a shortcoming of the GGB method; rather,
it is a reflection of the difficulty of the problem and the fact that the assumed sparsity level is
in this case high relative to n. Even when B ∗ = 0 (i.e., Σ∗ is a diagonal matrix), we still must
estimate p free parameters on the basis of only n observations. More generally, we know that the
rate of convergence in Theorem 1 cannot, up to log factors, be improved by any other method. To
see this, note that in the special case that G is a path graph, Bien et al. (2015a) prove (in Theorem
7 of their paper) that the optimal rate is B ∗ p/n over the class of B ∗ -banded matrices (assuming
√
B ∗ ≤ n). Since |gB ∗ | ∼ B ∗ p for a path graph, we see that the GGB estimator is within a log
factor of the optimal rate in Frobenius norm.
If we are interested in estimating the correlation matrix R∗ = diag(Σ∗ )−1/2 Σ∗ diag(Σ∗ )−1/2
instead of the covariance matrix, then we can apply our estimator to the sample correlation matrix
R̂ = diag(S)−1/2 · S · diag(S)−1/2 . Because the diagonal entries of the correlation matrix do not
need to be estimated, we can show (under one additional assumption) that the rate of convergence
has |gB ∗ | in place of p + |gB ∗ |.
Assumption 3 (bounded non-zero variances). Suppose that there exists a constant κ > 0 (i.e.,
not depending on n or p) such that minj Σ∗jj ≥ κ.
Theorem 2. Suppose Σ∗ is B ∗ -banded with respect to the graph G (in the sense of Definition 1)
and Assumptions
1, 2, and 3 hold. Then there exist positive constants x and c for which taking
q
λ=x

log(max{p,n})
,
n

wb = |gb |1/2 , δ = −∞, and M = diam(G),
kΣ̂λ (R̂) − R∗ k2F .

|gB ∗ | log(max{p, n})
n

with probability at least 1 − c/ max{p, n}.
When G has a small number of edges (e.g., with all but a small number of nodes having 0
degree), then we can have |gB ∗ |  n < p, in which case the above result establishes the consistency
of the graph-guided estimator even when p > n.

8

2.2.2

Bandwidth recovery

The remaining properties of this section can be established under a condition on the minimum
signal strength (analogous to the “beta-min” conditions arising in the regression context).
Assumption 4 (minimum signal strength). The nonzero elements of Σ∗ must be sufficiently large:
min |Σ∗jk | > 2λ.

jk∈gB ∗

q

log(max{p,n})
One can think of λ as on the order of
and measuring the “noise level” of the
n
data, in particular kS − Σ∗ k∞ . Thus, this assumption ensures that the nonzeros stand out from the
noise. Under this additional assumption, our estimator’s bandwidth, B(Σ̂λ ), matches the unknown
bandwidth B ∗ with high probability.

Theorem 3. Adding Assumption 4 to the conditions of Theorem 1, there exists a constant c > 0
such that B(Σ̂λ ) = B ∗ with probability at least 1 − c/ max{p, n}.
Proof. See Appendix C.
2.2.3

Positive semidefinite constraint

The results presented above are for δ = −∞. The next result provides an additional condition
under which the previous results hold for the δ = 0 case (i.e., the situation in which the estimator
is forced to be positive semidefinite).
Suppose all nonzero elements of Σ∗ are within a factor of τ of each other:
maxjk∈gB∗ |Σ∗jk |
≤ τ.
minjk∈gB∗ |Σ∗jk |
Assumption 5 (minimum eigenvalue). The minimum eigenvalue of Σ∗ must be sufficiently large:
p
λmin (Σ∗ ) > 4xτ log(p)/n max |gjB ∗ |.
j

The quantity maxj |gjB ∗ | is the maximal degree of the true covariance graph (defined by the
sparsity pattern of Σ∗ ).
Theorem 4. Adding Assumptions 4 and 5 to the conditions of Theorem 1, Theorems 1 and 3 hold
for δ = 0.

2.3

Discussion & the small world phenomenon

The theorems above assume that M = diam(G), but the results in fact apply more generally so
long as M ≥ B ∗ . While in general we may not know B ∗ , there are some applications in which
one may know enough about the mechanism of the graph-based dependence to provide an upper
bound on B ∗ . In such situations, choosing an M that is much smaller than diam(G) may lead to
substantial computational improvements (both in terms of timing and memory).
9

Examining the Frobenius-norm bound above, we see that the structure of the seed graph G plays
an important role in the rates of convergence. While it is tempting to think of B ∗ as controlling
the sparsity level (and therefore the difficulty of the problem), it is important to observe that this
quantity is in fact a function of both B ∗ and the graph G. For example, in the Frobenius-norm
bound, the quantity |gB ∗ | is the number of pairs of variables that are within B ∗ hops of each other in
the graph G. In conventional banded estimation, where G is a path seed graph, the dependence on
B ∗ is linear: |gB ∗ | ∼ B ∗ p. More generally, when G is a d-dimensional lattice graph |gB ∗ | ∼ (B ∗ )d p
(when B ∗ is small enough to avoid edge effects). However, intuition provided by lattice graphs does
not extend well to more
 general graphs. For a particularly extreme example, when G is a complete
p
graph, then |g1 | = 2 . Many large graphs in nature have small diameter (Watts & Strogatz 1998),
which means that |gB ∗ | may be quite large even when B ∗ is small. The idea of “six degrees of
separation” (Travers & Milgram 1969) means that when using a social network as a seed graph,
|g6 | ∼ p2 (in fact, recent work suggests six might even be an overestimate, Backstrom et al. 2011).
When diam(G) is very small, Definition 1 is restricted to a very limited class of sparsity
patterns. In the next section, we consider a much larger class of sparsity patterns, using the notion
of local bandwidths introduced in Definition 2.

3

Graph-guided banding with local bandwidths

Definition 2 describes a different kind of graph-guided banding, in which each variable j has an
associated bandwidth Bj∗ , which can be thought of as a neighborhood size. Consider sets of the
form
gjb = {k : 1 ≤ dG (j, k) ≤ b}, for 1 ≤ b ≤ Mj , 1 ≤ j ≤ p,
which is the set of variables within b hops of j in the seed graph G. We can describe a sparsity
pattern of the type in Definition 2 as
i h
i
[ h
S∗ =
{j} × gjBj∗ ∪ gjBj∗ × {j} .
1≤j≤p

To see this, note that jk ∈ S ∗ is equivalent to dG (j, k) ≤ max{Bj∗ , Bk∗ }. We desire an estimator
whose sparsity patterns are unions of the {gjb }, so we use the latent overlapping group lasso
penalty (Jacob et al. 2009), which is designed for such situations. We define the graph-guided
banding estimator with local bandwidths (GGB-local) as


1
2
Σ̃λ = arg min
kS − ΣkF + λP̃ (Σ; G) ,
(4)
ΣδIp 2
where
P̃ (Σ; G) =

min
{V(jb) ∈Rp×p }


Mj
p X
X


wjb kV(jb) kF

j=1 b=1

s.t.

ΣgjMj =

Mj
X
b=1

Vg(jb)
jb




.



The only difference between this estimator and the GGB-global estimator of (2) is in the choice
of groups. For a fixed j, the groups are nested, gj1 ⊆ · · · ⊆ gjMj . When Mj is the diameter of
the connected component that contains j, then gjMj represents all variables that are connected to
variable j. For j 6= k, we can have gjb ∩gkb 6= ∅ without one group contained completely in the other.
The fact that this group structure is not simply hierarchically nested, as in the global bandwidth
estimator, complicates both the computation and the theoretical analysis of the estimator.
10

3.1

Computation

By Lemma 1, we can solve (4) by Algorithm 1. The challenge of evaluating the GGB-local estimator
lies therefore in efficiently evaluating the proximal operator ProxλP̃ . Viewing this problem as an
optimization problem over the latent variables {V(jb) } suggests a simple blockwise coordinate
descent approach, which we adopt here.

3.2

Theory

In Theorem 1, we saw that the Frobenius-norm convergence depends on the number of nonzero
off-diagonal elements, |gB∗ |, a quantity that we observed in Section 2.3 could be quite large—
as large as O(p2 )—even when B ∗ is very small. In other words, for such “small-world” seed
graphs G, the assumption of graph-guided bandedness with a small global bandwidth does not
necessarily correspond to Σ∗ being a sparse matrix. The notion of graph-guided bandedness with
local bandwidths, as in Definition 2, provides finer control, allowing us to describe sparse Σ∗ even
with “small-world” seed graphs G. The purpose of the next theorem is to show that, under such
a notion of sparsity, the GGB-local estimator can attain Frobenius-norm consistency when Σ∗ is
sparse.
Theorem 5. Suppose Σ∗ is (B1∗ , . . . , Bp∗ )-banded with respect to a graph G (in the sense of Defq
inition 2) and Assumptions 1-2 hold. Then with λ = x log(max{p,n})
(for x sufficiently large),
n
wjb = |gjb |1/2 , δ = −∞, and Mj = maxk dG (j, k),
r
log(max{p, n}) p log(max{p, n})
∗ 2
∗
kΣ̃λ − Σ kF . |S |
+
n
n
with probability at least 1 − c/ max{p, n}.
Proof. See Appendix D.
The two terms on the right-hand side of Theorem 5 represent estimation of the
q off-diagonal and
∗
∗
∗
, then
diagonal elements of Σ . If Σ is sparse enough—in particular that |S | ≤ p log(max{p,n})
n
∗
estimating all of Σ is as difficult as estimating its diagonal elements alone (i.e., just the p variances).
A simple example of this situation is when Σ∗ is (B1∗ , . . . , Bp∗ )-banded with respect to a graph G
q
and Bj∗ > 0 for at most log(max{p,n})
variables j. Another example is when all Bj∗ ≤ 1 and the
n
q
sum of degrees of all j with Bj∗ = 1 is bounded by p log(max{p,n})
.
n
Theorem 5 suggests that consistency in Frobenius norm might not be achieved in the p > n
regime even in the sparsest situation of a diagonal matrix. As noted earlier, this is not a shortcoming
of our method but rather a reflection of the difficulty of the task, which requires estimating p
independent parameters even in the sparsest situation of a diagonal matrix.
As noted previously, if instead we were interested in estimating the correlation matrix, we can
attain a rate as in Theorem 5 but without the second term by applying our estimator to the sample
correlation matrix.
Theorem 6. Suppose Σ∗ is (B1∗ , . . . , Bp∗ )-banded with respect to a graph G (in the sense of Definition 2) and Assumptions 1, 2, and 3 hold. Then there exists a constant x for which taking
11

λ=x

q

log(max{p,n})
,
n

wjb = |gjb |1/2 , δ = −∞, and Mj = maxk dG (j, k),
r
kΣ̃λ (R̂) −

R∗ k2F

∗

. |S |

log(max{p, n})
n

with probability at least 1 − c/ max{p, n}.
Consistency follows even when p  n as long as the right-hand side approaches 0.

4
4.1

Empirical study
Simulations

All simulations presented in this paper were performed using the simulator package (Bien 2016b,a)
in R (R Core Team 2017). We investigate the advantages of exploiting a known graph structure
(both in the global and local bandwidths settings), and we also study our two methods’ robustness
to graph misspecification. We compare the two GGB methods to each other and to two other
methods, which represent the “patternless sparsity” and traditional banding approaches.
Perhaps the simplest and most commonly used sparse covariance estimator is to simply softthreshold each off-diagonal element of the sample covariance matrix:
soft

Σ̂jk = sign(Sjk ) [|Sjk | − λ]+ .
This method (Rothman et al. 2009), along with its close relatives (Rothman 2012, Liu et al. 2013,
Xue et al. 2012), do not make use of the graph G and corresponds to (3) with Ω being the `1 norm
and δ = −∞. Thus, comparisons to this method allow us to understand what is gained and what
is lost by using the graph G versus seeking “patternless sparsity.”
We also wish to compare our new estimators to a more traditional (i.e., non-graph-guided)
banding approach. Graph-guided banding is not the only way one might incorporate the graph
G. A different approach would be to find an ordering of the nodes based on G and then apply a
traditional banding estimator with respect to this ordering. To this end, we seek an ordering of the
variables that travels along G as much as it can, making jumps to another part of G as sparingly
as possible. We express this as a traveling salesman problem (TSP), which we then solve using an
off-the-shelf TSP solver (see Appendix F for details and for images of the learned paths). We then
apply the convex banding method using this ordering (Bien et al. 2015a,b).
While each of these methods can be solved with a positive semidefinite constraint, we compare
their performances without such a constraint. Doing so allows us to observe the fraction of the
time each unconstrained estimator is positive semidefinite: Table 5 in Appendix G shows that the
banding estimators are observed in the study to always be positive semidefinite, whereas this is not
the case for the soft-thresholding estimator. It follows that the reported results for the banding
methods in fact also pertain to these estimators with δ = 0.
To generate a covariance matrix that is (B1∗ , . . . , Bp∗ ) banded with respect to a graph G, we take
(
Σ∗jk =

1
dG (j,k) 1

n
o
dG (j, k) ≤ max(Bj∗ , Bk∗ )

a

if j 6= k
if j = k,

12

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

Figure 4: Two-dimensional lattice graph; seed graph (left panel), global bandwidth (center panel), local
bandwidths (right panel)
●

●

●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●
●

●
●

●

●

●

●
●

●

●
●

●

●

●
●

●

●

●

● ●

●

●
●

●

●

●

●

●
●

●

●

●

●

●

●

●
●
●

●

●

●

●

●
●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●
●
●

●

●

●
●

●

●

●

●

●

●
●

●

●

●

●

●

●
●

●

●
●

●

●

●
●

●
●

●

●
●
●

●

●

●

●

●

●

●
●

●

●
●
●

●

●

●

●

●

●

●

●
●
●

●
●

●
●

●
●

●

●

●

●

●

●

●
●
●

●

●
●

●

●

●

●
●

●
●

●
●

●
●

●
●

●

●

●

●
●

●

●

●

●
●

●

●

●

●
●

●

●

●
●

●
●

●

●

●

●

●

●

●
●

●
●

●

●

●
●
●

●

●
●

●

●

●

●

●

●

●

●

●

●
●

●

●

●
●

●
●

●
●

●

●

●

●
●

●

●

●
●

●

●

●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●
● ●

●
●

●

●

●

●

●

●
●

●

●
●

●
●
●

●

●

●

●
●

●

●

●

●

●

●

●
● ●
● ●

●

● ●
●
●
●

●

●

●

●

●
●

●

●
●

●
●

●

●
●

●
●

●

●
●

●
●
●

●

●
●

●

●

●

●

●

●
●

●
●

●
●

●

●

●

●

●
●

●
●

●

●

●

●

● ●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

● ●

●

●

●

●

●

●

●
●

●

●

●

●

●

● ●
● ●
●
●
● ●

●

●
●

●
●

●

●
●

●
●

●
●

●

●

●
●
●

●

●

●
●
●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

● ● ●
●
● ●
●
●
●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

●

●
●

●

●

●

●
●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●
●
●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●

●

●
●

●
●

●

●
●

●
●
●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●

● ●

●

●

●
●

●

●

●

●
●

●

●

●

●
●

●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●
●
●

●
●

●

●

●

●

●

●

●

●

●
●

●
●

● ●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●
●

●

●

●
●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●
●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●
●

●

●

●

●
●

●

● ●

●

●

●

●

●

●
●

●

●
●

●
●

●

●

●

●

Figure 5: Scale-free graph; seed graph (left panel), local bandwidths (right panel)
where a is chosen to ensure that the minimum eigenvalue of Σ is at least σ 2 (we take σ = 0.01
throughout).
Initially, we compare the methods under four GGB structures: two-dimensional lattice with
global bandwidth, two-dimensional lattice with local bandwidths, scale-free graph with global bandwidth, and scale-free graph with local bandwidths. In the first two scenarios, we take G to be a
two-dimensional square lattice having 20 variables per side (so p = 400). The global bandwidth is
taken to be B ∗ = 4, and the local bandwidths are generated according to


0 with probability 0.90
∗
Bj = 1 with probability 0.06


4 with probability 0.04
(see Figure 4). The resulting Σ∗ has roughly 10% of elements nonzero in the global case and about
1% in the local case. In the third and fourth scenarios, we take G to be a scale free network (as
generated by the function sample pa in the R package igraph, Csardi & Nepusz 2006) with p = 400
nodes. The same two choices of global and local bandwidths lead to Σ∗ having, respectively, about
75% and 10% nonzeros (see Figure 5). The marked difference in sparsity levels is a manifestation
of the differing neighborhood structures between these two graph types.

13

Table 1: A comparison of mean Frobenius error at CV-best (averaged over 200 replicates).
2-d lattice (20 by 20), b = 4
2-d lattice (20 by 20) w/ var-bw
Scale free (p = 400), b = 4
Scale free (p = 400) w/ var-bw

GGB-global
7.605 (0.015)
4.244 (0.002)
6.740 (0.009)
3.273 (0.002)

GGB-local
8.381 (0.031)
3.170 (0.008)
6.543 (0.015)
2.867 (0.007)

Soft thresholding
11.634 (0.014)
4.347 (0.007)
7.655 (0.004)
3.552 (0.005)

TSP + hierband
11.562 (0.031)
4.779 (0.003)
7.727 (0.009)
3.661 (0.002)

Table 2: A comparison of mean operator error at CV-best (averaged over 200 replicates).
2-d lattice (20 by 20), b = 4
2-d lattice (20 by 20) w/ var-bw
Scale free (p = 400), b = 4
Scale free (p = 400) w/ var-bw

4.1.1

GGB-global
2.185 (0.007)
0.989 (0.003)
4.350 (0.010)
1.160 (0.001)

GGB-local
2.428 (0.015)
0.756 (0.005)
4.110 (0.016)
1.008 (0.005)

Soft thresholding
3.224 (0.010)
0.997 (0.004)
4.695 (0.004)
1.212 (0.002)

TSP + hierband
3.149 (0.016)
1.024 (0.003)
4.709 (0.009)
1.238 (0.001)

Lattice versus scale-free graphs and global versus local bandwidths

We generate samples of size n = 300 under the four scenarios (and repeat 200 times). Tables 1 and
2 present Frobenius-norm error and operator-norm errors, kΣ̂−Σ∗ kF and kΣ̂−Σ∗ kop , respectively;
table 3 presents the entropy error, − log det(Σ̂Σ∗−1 ) + trace(Σ̂Σ∗−1 ) − p. Each method has a single
tuning parameter, and we report the performance of the method with the value of the tuning
parameter chosen by 5-fold cross validation. In particular, we randomly
partition the training data
P
into 5 folds and, for an estimator Σ̂λ , choose λ̂CV-best = arg minλ 51 5k=1 kS(k) −Σ̂λ (S(−k) )kF , where
S(k) is the sample covariance matrix in which only the kth fold is included and S(−k) is the sample
covariance matrix in which the kth fold is omitted. In the 2-d lattice model, GGB-global does best
when the true model has global bandwidth, as we would expect it to. Interestingly, GGB-local does
best not just when the true model has local bandwidths but also in the scale-free model with a global
bandwidth. This is a reminder of the important role that the seed graph’s structure plays in the
statistical properties of the estimators. Finally, Figure 6 shows ROC curves for the four methods.
The results are consistent with the results from the tables. In particular, GGB-global does best in
the first scenario, and GGB-local does best in the other scenarios. Interestingly, TSP + hiernet
does not appear to be an effective approach. Unlike the graph-guided approaches, TSP + hiernet
does worse than soft-thresholding in nearly every scenario. For example, in Figure 6, we see that
TSP + hiernet’s ability to identify nonzero elements of Σ∗ is no better than random guessing. This
poor performance is to be expected in this scenario since such a network cannot be well captured
by an ordering of the nodes. For the two-dimensional lattice graph scenarios, we find that such
an approach is at least partially successful, which makes sense since the ordering captures some
aspects of the underlying seed graph.

Table 3: A comparison of mean entropy error at CV-best (averaged over 200 replicates).
2-d lattice (20 by 20), b = 4
2-d lattice (20 by 20) w/ var-bw
Scale free (p = 400), b = 4
Scale free (p = 400) w/ var-bw

GGB-global
7719 ( 24)
25823 ( 55)
72606 (146)
52047 ( 70)

GGB-local
7819 ( 69)
18088 (119)
68087 (191)
44188 (186)

14

Soft thresholding
10234 ( 34)
25358 ( 89)
68654 ( 62)
53332 ( 86)

TSP + hierband
11870 ( 73)
27188 ( 62)
68518 (119)
54304 ( 76)

2−d lattice (20 by 20), b = 4

2−d lattice (20 by 20) w/ var−bw

1.00

0.75

0.50

0.25

Specificity

Method
0.00

GGB − global
Scale free (p = 400), b = 4

GGB − local

Scale free (p = 400) w/ var−bw

Soft thresholding

1.00

TSP + hierband

0.75

0.50

0.25

0.00
0.00

0.25

0.50

0.75

1.00 0.00

0.25

0.50

0.75

1.00

Sensitivity

Figure 6: ROC curves (formed by varying the tuning parameter for each method), show the sensitivity and specificity of identifying the nonzero elements of Σ∗ . For visibility, only five of the 200
replicates are shown. (Plot made using R package ggplot2, Wickham 2009.)

15

4.1.2

Effect of graph misspecification

The previous simulation studies establish the rather obvious fact that when one knows the true
graph G with respect to which Σ∗ is banded, one can better estimate Σ∗ with an estimator that
uses G. In practice, however, one might typically use a graph G0 that is (hopefully) close to G but
not exactly identical. This section seeks to understand the effect of graph misspecification on the
effectiveness of graph-guided banding. We simulate data as above but imagine that we only have
access to a corrupted version of the underlying graph G.
Let π ∈ [0, 1] control the degree of corruption to G. Of course there are many ways to corrupt a
graph, but we choose a scheme that will preserve the properties of the graph and the sparsity level
of the covariance matrix (so as not to confound graph distortion with other factors such as sparsity
level, number of triangles, etc.). In particular, we take G to be an Erdos-Renyi graph with p = 400
vertices and 760 edges (this is chosen to match the number of edges in the previous examples). We
take Σ∗ to have global graph-guided bandwidth B ∗ = 4 with respect to this graph G. For each
edge in G, with probability π we remove it and then connect two currently unconnected nodes
chosen at random. We choose G to be an Erdos-Renyi graph so that the distorted graph G0 is
itself Erdos-Renyi with the same number of edges. By keeping all graph properties fixed, we can
be confident that as we vary π, any differences in performance of the methods will be attributable
to G and G0 having different edges.
Figure 7 shows the performance of the three methods as the probability of reassigning an edge
is varied from 0 (no distortion of G) to 1 (all edges are distorted). As expected, we see both
GGB methods deteriorating as the distortion level of the seed graph G0 is increased. When π = 0
(i.e., G = G0 ) GGB-global is best (which makes sense since Σ∗ has a global bandwidth). When
approximately π = 0.25 of the edges are distorted, the best obtainable operator-norm error by
the GGB methods becomes worse than that of soft-thresholding. Soft-thresholding is, of course,
unaffected by graph distortion since it does not use G0 at all. Surprisingly, in Frobenius-norm error
the two GGB methods remain better than soft-thresholding even at very high levels of distortion.
A deeper examination of this simulation reveals that the GGB methods are more biased but less
variable than the soft-thresholding estimates, leading to a smaller Frobenius-norm error. We do
expect this behavior to be problem-specific and in general would expect GGB methods with highly
distorted graphs to perform worse than soft-thresholding.

4.2

Application to images of handwritten digits

In this section, we apply our estimators to a dataset of (16 × 16)-pixel images of handwritten digits
(Le Cun et al. 1990, Hastie et al. 2009)1 . There are between (roughly) 700 and 1500 images per
digit type. For each digit type, we take a random n = 90 images to form the “training set” sample
covariance matrix Strain and use the remaining data to form the “test set” sample covariance matrix
Stest (our results are based on 20 replications of this process). A natural seed graph in this situation
is the two-dimensional lattice graph formed by connecting adjacent pixels. We compare our two
GGB estimators to the soft-thresholding estimator. We fit each covariance estimator to the training
data and compute its best attainable (over choice of tuning parameter) distance in Frobenius norm
to Stest ,
min kΣ̂λ (Strain ) − Stest k2F .
λ≥0

1

available at http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/zip.digits/

16

6
●

●

●

15

●

●

●
●
●

●

●

●
●

●
●

●
●
●

5

●
●

●
●
●

●
●
●

●
●

●
●
●

●

●
●

●

3

4

●
●
●

2

Best operator error

10

●
●

●
●

1

5

Best Frobenius error

●
●

0

●
●

0.0

GGB − global
GGB − local
Soft thresholding
0.2

0

●

0.4

0.6

0.8

1.0

0.0

Probability of distortion

0.2

0.4

0.6

0.8

1.0

Probability of distortion

Figure 7: The effect of graph misspecification on the performance of the global and local graph-guided banding
methods in terms of Frobenius and operator norms.

Table 4: A comparison of best attainable Frobenius-norm test error (averaged over 20 replicates). Standard
errors are reported in parentheses. For all but two of the digits, a GGB method has lower test error than
soft-thresholding in all 20 out of 20 replications.
Digit
Digit
Digit
Digit
Digit
Digit
Digit
Digit
Digit
Digit

0
1
2
3
4
5
6
7
8
9

GGB-global
10.8 (0.1)
2.6 (0.1)
10.9 (0.2)
8.9 (0.1)
9.4 (0.1)
10.7 (0.2)
9.2 (0.2)
7.3 (0.1)
9.3 (0.1)
7.8 (0.2)

GGB-local
10.8 (0.1)
2.6 (0.1)
10.9 (0.2)
8.9 (0.1)
9.5 (0.1)
10.8 (0.2)
9.1 (0.2)
7.4 (0.1)
9.4 (0.1)
7.9 (0.2)

Soft thresholding
11.0 (0.1)
2.6 (0.1)
11.3 (0.1)
9.3 (0.1)
9.9 (0.1)
11.2 (0.2)
9.2 (0.2)
7.5 (0.1)
9.8 (0.1)
8.0 (0.2)

Table 4 reports this quantity, averaged over the 20 random train-test splits of the dataset, for
each digit. GGB-local has a lower test error than soft-thresholding in all 20 replicates for the
digits 2-9; GGB-global has a lower test error than soft-thresholding in all 20 replicates for all digits
other than 0, 1, and 6. In fact, there is no digit for which soft-thresholding wins in a majority
of replicates. This indicates that exploiting the neighborhood structure inherent to image data is
useful for covariance estimates in this context.

5

Discussion

This paper generalizes the notion of banded covariance estimation to situations in which there is
a known underlying graph structure to the variables being measured. The traditional notion of a
banded matrix corresponds to the special case where this graph is a path encoding the ordering
of the variables. We propose two new “graph-guided” estimators, based on convex optimization
17

problems: GGB-global enforces a single bandwidth across all variables while GGB-local allows each
variable to have its own bandwidth. There are many interesting avenues of future research. In this
paper, we assume the seed graph is given, and the challenge is in determining the bandwidth(s).
One could imagine situations in which the seed graph itself is unavailable (but in which we do
believe that such a graph exists). In such a setting, can one identify the smallest seed graph G for
which the covariance matrix’s sparsity pattern is banded with respect to G? This generalizes the
problem of discovering an ordering for which one can apply a traditional known-ordering estimator
considered by Wagaman & Levina (2009), Zhu & Liu (2009). A simple heuristic strategy might
involve forming a minimum spanning tree based on the complete graph with edge weights given
by the sample covariance matrix. The graph-guided banding estimators would then have sparsity
patterns with this minimum spanning tree seed graph serving as backbone. Another future avenue
of research is in applying the graph-guided structure to the inverse covariance matrix. Finally,
we have seen (e.g., in Theorem 1 and in the empirical study) that the properties of the seed
graph G influence the success of the GGB estimators. More generally understanding how the
graph properties of G affect the difficulty of the estimation problem itself would be a particularly
interesting direction of future study.

A

Proof of Lemma 1 of main paper

Lemma 2. Consider the problem

Σ̂ = arg min

ΣδIp


1
2
kS − ΣkF + λΩ(Σ) ,
2

where Ω is a closed, proper convex function, and let


1
2
ProxλΩ (M) = arg min
kM − ΣkF + λΩ(Σ)
Σ
2
denote the proximal operator of λΩ(·) evaluated at a matrix M. Algorithm 1 converges to Σ̂.
Proof. Problem (3) is equivalent to

1
min
kS − Σk2F + λΩ(Σ̃)
2
Σ,Σ̃


s.t. Σ = Σ̃,

Σ  δIp ,

which has Lagrangian
1
L(Σ, Σ̃; B, C) = kS − Σk2F + λΩ(Σ̃) + hΣ − Σ̃, Bi + hδIp − Σ, Ci.
2
The dual function is given by minimizing over Σ and Σ̃:


n
o
1
inf L(Σ, Σ̃; B, C) = δhIp , Ci + inf
kS − Σk2F + hΣ, B − Ci + inf λΩ(Σ̃) − hΣ̃, Bi
Σ
2
Σ,Σ̃
Σ̃
1
1
= δhIp , Ci − kB − C − Sk2F + kSk2F − [λΩ]∗ (B)
2
2

18

where [λΩ]∗ is the convex conjugate of λΩ. Also,
Σ̂ = S − B̂ + Ĉ.
The dual optimal (B̂, Ĉ) is given by solving
1
min kB − C − Sk2F + [λΩ]∗ (B) − δhIp , Ci s.t.
B,C 2

C  0,

which is a convex problem. Furthermore, by Tseng (2001), it can be solved by alternately optimizing
over B (holding C fixed) and over C (holding B fixed). The solution to the optimization over B is
the proximal operator of [λΩ]∗ evaluated at C + S. The Moreau decomposition (see, e.g., Parikh
& Boyd 2014) for the closed, proper convex function λΩ gives us
B̂(C) = Prox[λΩ]∗ (S + C) = S + C − ProxλΩ (S + C).
Optimizing over C while holding B fixed requires solving
1
Ĉ(B) = min kB − C − Sk2F − δhIp , Ci
C 2

s.t. C  0,

which in terms of the eigenvalue decomposition of B − S = UΛUT becomes
1
Ĉ(B) = min kC − UΛUT k2F − δhIp , Ci s.t. C  0
C 2
or Ĉ(B) = UD̂(B)UT where
1
D̂(B) = min kD − Λk2F − δhIp , Di
D 2

s.t. D  0,

which is solved by
D̂(B) = diag([Λii + δ]+ ).

B

Proof of Theorem 1 of main paper (Frobenius error)

Proof. Since δ = −∞ and M = diam(G), the estimator is more simply


1
2
Σ̂ := Σ̂λ = arg min
kS − ΣkF + λP (Σ; G)
Σ
2
where
P (Σ; G) =

min
{V(b) ∈Rp×p }

(M
X

wb kV(b) kF

b=1

s.t.

Σ− =

M
X

(5)

)
Vg(b)
b

.

b=1

In what follows, whenever we write a sum with subscript b, it is understood to run along b =
1, . . . , M ; this also holds for the {V(b) } notation, which is short for {V(b) }M
b=1 .

19

We begin by observing that (5) is a strongly convex problem in Σ, and thus there is a unique
minimizer Σ̂, although the minimizing set of V̂(b) ’s may not be unique. Also, since diag(Σ) does
not appear in the penalty, we see that diag(Σ̂) = diag(S). Thus, we can write (5) as


1 −
−
2
kS − ΣkF + λP (Σ; G) .
(6)
Σ̂ = arg min
Σ
2
We can write (6) in terms of the V(b) ’s:
(
X
X
1
{V̂(b) } ∈ arg min
kS −
V(b) k2F + λ
wb kV(b) kF
2
V(b)
b

)
s.t.

(b)
Vg c
b

= 0 ∀b .

(7)

b

Note that we have replaced S− by S (for notational convenience) since gb does not include diagonal
elements. The lemma below provides necessary and sufficient conditions for {V̂(b) } to be a solution
to (7).
P
(b)
Lemma 3. Writing Σ̂ = b V̂(b) , {V̂(b) } is a solution to (7) iff. for each b ∈ {1, . . . , M }, V̂gc = 0
b
and
(b)
λwb V̂gb
[S − Σ̂]gb =
if V̂(b) 6= 0
kV̂(b) kF
and k[S − Σ̂]gb kF ≤ λwb otherwise.
Proof. Problem (7) is equivalent to a non-overlapping group lasso problem, and this lemma provides
the KKT conditions of this problem, which are well-known.
Corollary 1. Let {V̂(b) } be a solution to (7). Then,
k[S − Σ̂]gb kF < λwb =⇒ V̂(b) = 0.
Suppose supp(Σ∗ ) = gB ∗ . If we can construct a solution {V̂(b) } to (7) that has
max wb−1 k[S − Σ̂]gb kF < λ,

b>B ∗

then, by Corollary 1, we know that V̂(b) = 0 for all b > B ∗ or
{b : V̂(b) 6= 0} ⊆ {1, . . . , B ∗ }.
In fact, the corollary implies that this holds for all solutions, not just for our constructed {V̂(b) }.
Consider the restricted problem
(
)
B∗
B∗
X
X
1
(b)
(b) 2
(b)
(b) B ∗
∗
kS −
V kF + λ
wb kV kF s.t. Vgc = 0 for 1 ≤ b ≤ B ,
{V̄ }b=1 ∈ arg min
b
2
V(b)
b=1

b=1

PB ∗

(b)
and write Σ̄ =
b=1 V̄ . This restricted problem is of an identical form to (7) but with M
replaced by B ∗ . Thus, Lemma 3 applies if we replace b = 1, . . . , M with b = 1, . . . , B ∗ : For each
b ≤ B∗,
(b)
λwb V̄gb
if V̄(b) 6= 0
[S − Σ̄]gb =
kV̄(b) kF

20

and k[S − Σ̄]gb kF ≤ λwb otherwise.
P
(b) =
Now, for b > B ∗ , let V̄(b) = 0. We show that {V̄(b) }M
is a solution to (7). Clearly, M
b=1 V̄
b=1
PB ∗ (b)
= Σ̄, so by Lemma 3 it simply remains to show that
b=1 V̄
k[S − Σ̄]gb kF ≤ λwb for b > B ∗ .
For b > B ∗ :
k(S − Σ̄)gb k2F = k(S − Σ̄)gB∗ k2F + kSs(B∗ +1):b k2F
2
2
∗
≤ λ2 wB
∗ + k(S − Σ )s(B ∗ +1):b kF
2
∗ 2
≤ λ2 wB
∗ + |s(B ∗ +1):b | · kS − Σ k∞ .

Now, on the event
Aλ = {kS − Σ∗ k∞ < λ},

(8)

we have
2
2 2
k(S − Σ̄)gb k2F < λ2 (wB
∗ + |s(B ∗ +1):b |) = λ wb ,

using that wb2 = |gb | = |gB ∗ | + |s(B ∗ +1):b |. This establishes, by Lemma 3, that {V̄(b) } is a solution
to (7) and, by Corollary 1 that {b : V̂(b) 6= 0} ⊆ {1, . . . , B ∗ } for all solutions {V̂(b) } to (7). Given
the support constraints of each V̂(b) , it follows that the overall Σ̂ has bandwidth no larger than
B∗.
Proposition 1. If wb = |gb |1/2 , δ = −∞, and M = diam(G), then B(Σ̂λ ) ≤ B ∗ on Aλ .
Having established the above proposition about the bandwidth, it is straightforward to prove a
result about the estimation error. On Aλ ,
k(Σ̂ − Σ∗ )gB∗ kF ≤ k(S − Σ∗ )gB∗ kF + k(S − Σ̂)gB∗ kF ≤ kS − Σ∗ k∞ · |gB ∗ |1/2 + λwB ∗ < 2λwB ∗
and Σ̂gB∗ = 0. Thus,
2
kΣ̂ − Σ∗ k2F = kdiag(S) − diag(Σ∗ )k2F + k(Σ̂ − Σ∗ )gB∗ k2F < λ2 (p + 4wB
∗ ).

(9)

Proposition 2. If wb = |gb |1/2 , δ = −∞, and M = diam(G), then
kΣ̂ − Σ∗ k2F ≤ λ2 (p + 4|gB(Σ∗ ) |)
on Aλ .
Under Assumptions 1 and 2 of the main paper, there exists a constant c > 0 for which
P (Aλ ) ≥ 1 − c/ max{p, n}

(10)

p
if we take λ = x log max{p, n}/n with x sufficiently large. This follows directly from bound
(A.4.1) of Theorem A.1 of Bien et al. (2015c). Theorem 1 of the main paper follows.

21

C

Proof of Theorem 3 of main paper (bandwidth recovery)

We prove a result that is a bit more general than Theorem 3 of the main paper.
Theorem 1. Under the minimum signal condition
min kΣ∗gb kF /|gb |1/2 > 2λ,

b≤B ∗

and the conditions of Theorem 1 of the main paper, B(Σ̂) = B(Σ∗ ) with probability at least 1 −
c/ max{p, n}.
Proof. On the event Aλ defined in (8), for b ≤ B ∗ ,
k(Σ̂ − Σ∗ )gb kF ≤ k(S − Σ∗ )gb kF + k(S − Σ̂)gb kF
< λ|gb |1/2 + k(S − Σ̂)gb kF
≤ λ|gb |1/2 + λwb = 2λwb .
Thus,
max wb−1 k(Σ̂ − Σ∗ )gb kF < 2λ.

b≤B ∗

We already know that B(Σ̂) ≤ B ∗ , so it suffices to show that
min kΣ̂gb kF /wb > 0,

b≤B ∗

which follows since
kΣ̂gb kF /wb ≥ kΣ∗gb kF /wb − k(Σ̂ − Σ∗ )gb kF /wb > 2λ − 2λ = 0.
The result then follows from (10).
Corollary 2. Under the minimum signal condition
min |Σ∗jk | > 2λ,

jk∈gB ∗

and the conditions of Theorem 1 of the main paper, B(Σ̂) = B(Σ∗ ) with probability at least 1 −
c/ max{p, n}.
Proof. For any b ≤ B ∗ ,
kΣ∗gb k2F /|gb | =

X

[Σ∗jk ]2 /|gb | >

jk∈gb

X
jk∈gb

Thus, the conditions of Theorem 1 hold.

22

4λ2 /|gb | = 4λ2 .

D

Proof of Theorem 5 (Frobenius error of local bandwidths estimator)

Proof. Let Σ̃ := Σ̃λ be the minimizer of Equation (4) of the main paper. By definition,
1
1
kS − Σ̃k2F + λP̃ (Σ̃; G) ≤ kS − Σ∗ k2F + λP̃ (Σ∗ ; G).
2
2
Since
kS − Σ̃k2F = kS − Σ∗ k2F + kΣ̃ − Σ∗ k2F − 2hS − Σ∗ , Σ̃ − Σ∗ i,
we have
1
kΣ̃ − Σ∗ k2F + λP̃ (Σ̃; G) ≤ hS − Σ∗ , Σ̃ − Σ∗ i + λP̃ (Σ∗ ; G).
2

(11)

−1
By Lemma 3 of Obozinski et al. (2011), the dual norm is P̃ ∗ (B; G) = maxjb wjb
kBgjb kF , and so

hS − Σ∗ , Σ̃ − Σ∗ i ≤ λ2 p + P̃ ∗ (S − Σ∗ ; G)P̃ (Σ̃ − Σ∗ ; G).

(12)

Define the index sets I = {jb : 1 ≤ b ≤ Mj } and I ∗ = {jb : 1 ≤ b ≤ Bj∗ }.
Observe that with wjb = |gjb |1/2
−1
kBgjb kF ≤ kBk∞ .
max wjb
jb∈I

It follows that
hS − Σ∗ , Σ̃ − Σ∗ i < λ2 p + λP̃ (Σ̃ − Σ∗ ; G),
as long as kS − Σ∗ k∞ < λ, that is on the event Aλ , as defined in (8). Thus, combining this with
(11),
h
i
1
kΣ̃ − Σ∗ k2F + λP̃ (Σ̃; G) ≤ λ2 p + λP̃ (Σ̃ − Σ∗ ; G) + λP̃ (Σ∗ ; G).
2
holds on Aλ . Using the triangle inequality (P̃ is a norm),
h
i
1
kΣ̃ − Σ∗ k2F ≤ λ2 p + λ P̃ (Σ̃ − Σ∗ ; G) − P̃ (Σ̃; G) + P̃ (Σ∗ ; G) ≤ λ2 p + 2λP̃ (Σ∗ ; G).
2
By definition of P̃ ,
P̃ (Σ∗ ; G) ≤

X

wjb kV(jb) kF

(13)

jb∈I
(jb)

for any choice of {V(jb) } such that Vgc = 0 and Σ∗ =
(j 0 b0 )

Consider taking Vjk

P

jb V

jb

(jb) .

= 0, for all jk ∈ {1, . . . , p}2 and j 0 b0 ∈ I, unless
(jBj∗ )

(i) dG (j, k) ≤ min{Bj∗ , Bk∗ }, in which case we take Vjk
(kBk∗ )

(ii) Bj∗ < dG (j, k) ≤ Bk∗ , in which case we take Vjk
23

(jBj∗ )

= Vkj
(kBk∗ )

= Vkj

(kBk∗ )

= Vjk

= Σ∗jk .

(kBk∗ )

= Vkj

= 12 Σ∗jk .

(jBj∗ )

(iii) Bk∗ < dG (j, k) ≤ Bj∗ , in which case we take Vjk

(jBj∗ )

= Vkj

= Σ∗jk .

We verify that this choice of {V(jb) } satisfies the necessary constraints. Clearly, V(jb) = 0 for
b 6= Bj∗ . If jk 6∈ gjBj∗ , then dG (j, k) > Bj∗ , which excludes cases (i) and (ii), the only two cases in
(jBj∗ )

(jBj∗ )

(jB ∗ )

) could be nonzero. Thus, Vgc ∗j = 0.
jB
j
P
0 b0 )
∗
(j
We next verify that j 0 b0 ∈I V
= Σ . In case (i),

which Vjk

(and Vkj

(jBj∗ )

0 0

X

[V(j b ) ]jk = Vjk

(kBk∗ )

= Σ∗jk .

+ Vjk

j 0 b0 ∈I

In case (ii),
X

0 0

(kBk∗ )

0 0

(jBj∗ )

[V(j b ) ]jk = Vjk

= Σ∗jk .

j 0 b0 ∈I

In case (iii),
X

[V(j b ) ]jk = Vjk

= Σ∗jk .

j 0 b0 ∈I
(jBj∗ )

Now, we have that |Vjk
X

∗

| ≤ |Σ∗jk |, which implies that kV(jBj ) kF ≤ kΣ∗gjB∗ kF . Thus,
j

wjb kV(jb) kF =

X

wjBj∗ kV

(jBj∗ )

kF ≤

j

jb∈I

X

wjBj∗ kΣ∗gjB∗ kF ≤ kΣ∗ k∞

j

j

X

|gjBj∗ |.

j

Now, each jk ∈ S ∗ appears only in gjBj∗ or gkBk∗ (or potentially both). Thus, by (13),
2|S ∗ |.

∗

∗

It follows that P̃ (Σ ; G) ≤ 2kΣ k∞

|S ∗ |,

P

j

|gjBj∗ | ≤

and so (on Aλ )

kΣ̃ − Σ∗ k2F ≤ λ2 p + 4kΣ∗ k∞ |S ∗ |λ.
Assumptions 1 and 2 of the main paper are used (as in Appendix B) to show that Aλ holds with
the desired probability.

E

Proofs about GGB on sample correlation matrix

We begin by establishing a concentration inequality for the correlation matrix that will be used
in the proofs of both of these theorems. This result is similar to Lemma A.2 of Liu et al. (2017)
except that we assume sub-Gaussian rather than Gaussian data. It is also similar to Lemma D.5
of Sun et al. (2017) except that we assume that S has been centered by the sample mean.
Lemma 4. Suppose that Assumptions 1, 2, and 3 hold. There exists a constant x (i.e., a value
not depending on n or p) such that for all n and p for which log{p, n}/n < 1/x2 ,
p
P(kR̂ − R∗ k∞ > x log(max{p, n})/n) ≤ 16/ max{p, n}.

24

Proof. The following proof is based on the proof of Lemma A.2 of Liu et al. (2017). Certain
details have been fleshed out and some modifications have been made to accommodate the differing
assumptions.
Let Aj = {|Sjj − Σ∗jj | ≤ Σ∗jj /2}. Observing that Aj = {1/2 ≤ Sjj /Σ∗jj ≤ 3/2}, it follows that
Sjj Skk
≤ (3/2)2
Σ∗jj Σ∗kk


Sjj Skk 1/2
=⇒ 1/2 ≤
≤ 3/2
Σ∗jj Σ∗kk

Aj ∩ Ak =⇒ (1/2)2 ≤

1/2

1/2

1/2

=⇒ (1/2)(Σ∗jj Σ∗kk )jk ≤ (Sjj Skk )jk ≤ (3/2)(Σ∗jj Σ∗kk )jk
1/2

1/2

1/2

=⇒ |(Sjj Skk )jk − (Σ∗jj Σ∗kk )jk | ≤ (1/2)(Σ∗jj Σ∗kk )jk .

(14)

Now,
{|R̂jk − R∗jk | > t} = {|Sjk − Σ∗jk | > t(Sjj Skk )1/2 }
⊆ {|Sjk − Σ∗jk | > t(Sjj Skk )1/2 } ∩ (Aj ∩ Ak )

[
(Aj ∩ Ak )c .

By (14) and the triangle inequality,
{|Sjk − Σ∗jk | > t(Sjj Skk )1/2 } ∩ (Ai ∩ Aj ) ⊆ {|Sjk − Σ∗jk | > (t/2)(Σ∗jj Σ∗kk )1/2 },
so
{|R̂jk − R∗jk | > t} ⊆ {|Sjk − Σ∗jk | > (t/2)(Σ∗jj Σ∗kk )1/2 }

[
(Aj ∩ Ak )c .

Now,
{max |R̂jk − R∗jk | > t} ⊆
j6=k

i
[h
{|Sjk − Σ∗jk | > (t/2)(Σ∗jj Σ∗kk )1/2 } ∪ Acj ∪ Ack
j6=k

=

[

{|Sjk − Σ∗jk | > (t/2)(Σ∗jj Σ∗kk )1/2 } ∪

⊆

Acj

j

j6=k

[

[

{|Sjk − Σ∗jk | > (t/2)κ} ∪

[

Acj

j

j6=k

⊆ {max |Sjk −
j6=k

Σ∗jk |

> (t/2)κ} ∪ {max |Sjj − Σ∗jj | > (1/2)κ}
j

∗

⊆ {kS − Σ k∞ > (κ/2) min{t, 1}}.
Since R̂jj = 1 = R∗jj , maxj6=k |R̂jk − R∗jk | = kR̂ − R∗ k∞ and thus we have established that
P(kR̂ − R∗ k∞ > t) ≤ P(kS − Σ∗ k∞ > (κ/2) min{t, 1}).
By Lemma A.1 of Bien et al. (2015c) (with a typo corrected), there exist constants c1 , c2 such
that
2
P(kS − Σ∗ k∞ > ) ≤ 2p2 e−c2 n(/κ̄) + 8pe−c1 n/κ̄
for any 0 <  < 2κ̄. This constraint on  implies that (/κ̄)2 < 2(/κ̄) and so
2

2

2p2 e−c2 n(/κ̄) + 8pe−c1 n/κ̄ < 2p2 e−c2 n(/κ̄) + 8pe−c1 n(/κ̄)
25

2 /2

.

Taking c3 = min{c2 , c1 /2},
2

2

P(kS − Σ∗ k∞ > ) ≤ (2p2 + 8p)e−c3 n(/κ̄) ≤ 16p2 e−c3 n(/κ̄) .
Observing that  = (κ/2) min{t, 1} < κ̄ < 2κ̄, we have that
2 /(2κ̄)2

P(kR̂ − R∗ k∞ > t) ≤ 16p2 e−c3 nκ

min{t2 ,1}

.
p
We take t = x log(max{p, n})/n for x = 12(κ̄/κ)2 /c3 . By Assumption 3, κ and κ̄ are treated as
constants. If x2 log(max{p, n})/n then min{t2 , 1} = t2 and
2 /(2κ̄)2 x2

P(kR̂ − R∗ k∞ > t) ≤ 16p2 e−c3 κ

log(max{p,n})

2 /(2κ̄)2 x2

≤ 16 max{p, n}2−c3 κ

.

Thus, we just need that 2 − c3 κ2 /(2κ̄)2 x2 ≤ −1, which can be ensured by choosing
x = 12(κ̄/κ)2 /c3 .

We define the event
Bλ = {kR̂ − R∗ k∞ < λ}.
Lemma 4 establishes that there exist constants x and c for which choosing λ = x
ensures that P(Bλ ) ≥ 1 − c/ max{n, p}.

E.1

(15)
q

log(max{p,n})
n

Proof of Theorem 2 (Frobenius error of global GGB on correlation matrix)

Let Σ̂λ (R̂) be the minimizer of Equation (2) of the main paper where S is replaced by R̂.
Proof of Theorem 2. The proof is nearly identical to that of Theorem 1 given in Section C. We
take Σ̂ to be Σ̂λ (R̂), we replace S with R̂, and we replace Aλ with Bλ . Instead of (9), we get
2
kR̂ − R∗ k2F = k(Σ̂ − Σ∗ )gB∗ k2F < 4λ2 wB
∗,

which leads to the stated result.

E.2

Proof of Theorem 6 (Frobenius error of local GGB on correlation matrix)

Let Σ̃λ (R̂) be the minimizer of Equation (4) of the main paper where S is replaced by R̂.
Proof of Theorem 6. The proof is nearly identical to that of Theorem 5 given in Section D. We
take Σ̃ to be Σ̃λ (R̂) and replace S with R̂. Equation 12 becomes instead
hR̂ − R∗ , Σ̃ − R∗ i ≤ P̃ ∗ (R̂ − R∗ ; G)P̃ (Σ̃ − R∗ ; G).
Following the same steps as in that proof, we note that on the set Bλ , defined in (15),
kΣ̃ − Σ∗ k2F ≤ 4kΣ∗ k∞ |S ∗ |λ.

26

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●
●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
● ●
●
●

●

●

●
●

●

● ●

●
●

●

●
● ●

● ●
● ●
●
●
● ●
● ●

●

●

●

●
●

●

●

●

●

● ●
●
●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●
●

●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●
●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●
●
●
●
●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●
●

●

●
●

●

●

●

● ●
●
●

●

●

●

●

●

●
●

● ●

●
●

●

●

●
●

●

●
●

●
●

●
●

●
●

●

●
●

●
●

●
●
●

● ●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ● ●

●

●

●
●

●
●

●

●

●

● ●

●
●

●

●

●

●
●

●
● ●

●
●

●

●
●
●

●

●
●
●

● ● ● ●
●
● ●

●

●

●

●

●

●

●

●
●
● ●

●
●

●

●

●

●

●

●

●

●

●
●
●

●
●
●

●
●

●

●
●

●

●
●
●

● ●
●
●

●

●

●

●
●

●

●
●

●
●

●
●

●

●

●

●

●

●

●

● ●
●

●

●

●

●

●

●
●

●
● ●
●
●
●

●

●

●
●

●
●

●

●
●

●

●

●

●

●●

●

●

●

●

●

●

● ●●
●

●

●

● ●

●

●

●

●

●

●

●

●
● ●

●
●

●

●

●

●

●

●

●

●

●

● ●
●

●

●
●

●
●

●
●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

● ●

●
●

●

●
●

●

●

●

●
●

●

●
●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●
● ●
● ● ●

●

●

●

●

●

●
●

●
●
●

●
●

●

●
●

● ●
●

●

●
●

●
● ●

●

●

●
●

●

●

●

● ●
●
●
●
● ●
● ● ●
● ●
●
●

●

Figure 8: TSP solutions used for the two graphs used in the simulation studies

F

Traveling salesman problem

Let C be a p × p matrix of costs, with Cjk being the cost of traveling from node j to node k. The
traveling salesman problem (TSP) seeks a loop v1 → v2 → . . . vp → v1 through all the nodes such
P
that p−1
`=1 Cv` v`+1 + Cvp v1 is minimal.
In our context, we have a graph G and wish to find an ordering of the nodes that is strongly
informed by G. In particular, we would like this path to travel along edges of G as much as possible
and when this is not possible and the path must “fly” between two nodes that are not connected
by a node, that it take the shortest flight possible. We therefore take
Cjk = dG (j, k) − 1.
We use the R package TSP (Hahsler & Hornik 2017), which provides an interface to the Concorde
solver (Applegate et al. 2017), and Figure 8 shows the resulting paths.

G

An empirical study of positive semidefiniteness

For each of the models considered in Section 4.1.1 of the main paper, we record the proportion
of the time each method is positive semidefinite. In particular, this average is over both the 200
simulation replicates and the values of the tuning parameter of each method. Table 5 shows that
throughout this study the banding methods (GGB and hierband) are always observed to give
positive semidefinite estimates. This is not the case for the soft-thresholding estimator.

H

A higher dimensional example

In Section 4.1.1, we study four models in which p = 400 and n = 300. The covariance matrix has
p(p + 1)/2 ≈ 80, 000 parameters. Thus, the number of parameters to be estimated greatly exceeds
27

Table 5: A comparison of propor PSD (averaged over a grid of tuning parameters and 200 replicates).
2-d lattice (20 by 20), b = 4
2-d lattice (20 by 20) w/ var-bw
Scale free (p = 400), b = 4
Scale free (p = 400) w/ var-bw

GGB
1.000
1.000
1.000
1.000

- global
(0.000)
(0.000)
(0.000)
(0.000)

GGB
1.000
1.000
1.000
1.000

- local
(0.000)
(0.000)
(0.000)
(0.000)

1.0
0.4

0.6

0.8

0.6

0.8
0.2

0.2
0.0

0.0
0.0

0.4

Specificity

0.8
0.6
0.4

GGB − global bandwidth
GGB − local bandwidth
Soft thresholding

0.2

Specificity

TSP + hierband
1.000 (0.000)
1.000 (0.000)
1.000 (0.000)
1.000 (0.000)

2−d lattice (20 by 20)

1.0

2−d lattice (50 by 50)

Soft thresholding
0.349 (0.008)
0.625 (0.008)
0.998 (0.001)
1.000 (0.000)

1.0

0.0

Sensitivity

0.2

0.4

0.6

0.8

1.0

Sensitivity

Figure 9: ROC curves for three methods when p = 400 and p = 2, 500.
the number of samples. Nonetheless, we can ask what happens if p is even larger than this. We
revisit the example of the 2-d lattice, increasing the side length from 20 to 50 so that p = 2, 500. In
this case the covariance matrix has p(p + 1)/2 ≈ 3, 000, 000 parameters. We focus in particular on
the GGB methods and soft thresholding. Figure 9 provides a side-by-side comparison of the ROC
curves (five replicates are shown per method) when one increases from p = 400 to p = 2, 500. The
curves in the two plots are different although this is difficult to discern by eye. This suggests that
we are effectively in the same regime in these two simulated scenarios.

References
Applegate, D., Bixby, R., Chvatal, V. & Cook, W. (2017), Concorde TSP Solver.
URL: http://www.math.uwaterloo.ca/tsp/concorde/
Backstrom, L., Boldi, P., Rosa, M., Ugander, J. & Vigna, S. (2011), ‘Four degrees of separation.
doi: arxiv1111. 4570v1’.
Banerjee, O., El Ghaoui, L. E. & d’Aspremont, A. (2008), ‘Model selection through sparse maximum
likelihood estimation for multivariate gaussian or binary data’, Journal of Machine Learning
Research 9, 485–516.
Bickel, P. J. & Levina, E. (2008a), ‘Regularized estimation of large covariance matrices’, Annals of
Statistics pp. 199–227.
28

Bickel, P. & Levina, E. (2008b), ‘Covariance regularization by thresholding’, Annals of Statistics
36, 2577–2604.
Bickel, P. & Lindner, M. (2012), ‘Approximating the inverse of banded matrices by banded matrices
with applications to probability and statistics’, Theory of Probability & Its Applications 56(1), 1–
20.
Bien, J. (2016a), simulator: An Engine for Running Simulations. R package version 0.2.0.
URL: https://CRAN.R-project.org/package=simulator
Bien, J. (2016b), ‘The simulator:
arXiv:1607.00021 .

An engine to streamline simulations’, arXiv preprint

Bien, J., Bunea, F. & Xiao, L. (2015a), ‘Convex banding of the covariance matrix’, Journal of the
American Statistical Association DOI: 10.1080/01621459.2015.1058265 .
Bien, J., Bunea, F. & Xiao, L. (2015b), hierband: Convex Banding of the Covariance Matrix. R
package version 1.0.
URL: https://CRAN.R-project.org/package=hierband
Bien, J., Bunea, F. & Xiao, L. (2015c), ‘Supplementary materials to “convex banding of the covariance matrix”’, Journal of the American Statistical Association .
Bien, J. & Tibshirani, R. (2011), ‘Sparse estimation of a covariance matrix’, Biometrika 98(4), 807–
820.
Cai, T., Liu, W. & Luo, X. (2011), ‘A constrained l1 minimization approach to sparse precision
matrix estimation’, Journal of the American Statistical Association 106(494), 594–607.
Cai, T. T. & Yuan, M. (2012), ‘Adaptive covariance matrix estimation through block thresholding’,
The Annals of Statistics 40(4), 2014–2042.
Cai, T. T. & Yuan, M. (2015), ‘Minimax and adaptive estimation of covariance operator for random
variables observed on a lattice graph’, Accepted to Journal of the American Statistical Association
.
Cai, T. T., Zhang, C.-H. & Zhou, H. H. (2010), ‘Optimal rates of convergence for covariance matrix
estimation’, The Annals of Statistics 38(4), 2118–2144.
Chandrasekaran, V., Parrilo, P. A. & Willsky, A. S. (2010), Latent variable graphical model selection
via convex optimization, in ‘Communication, Control, and Computing (Allerton), 2010 48th
Annual Allerton Conference on’, IEEE, pp. 1610–1613.
Csardi, G. & Nepusz, T. (2006), ‘The igraph software package for complex network research’,
InterJournal Complex Systems, 1695.
URL: http://igraph.org
Das, J. & Yu, H. (2012), ‘Hint: High-quality protein interactomes and their applications in understanding human disease’, BMC systems biology 6(1), 92.
Dempster, A. P. (1972), ‘Covariance selection’, Biometrics 28(1), 157–175.
29

Deng, X. & Yuan, M. (2009), ‘Large gaussian covariance matrix estimation with markov structures’,
Journal of Computational and Graphical Statistics 18(3).
Donoho, D. L., Gavish, M. & Johnstone, I. M. (2013), ‘Optimal shrinkage of eigenvalues in the
spiked covariance model’, arXiv preprint arXiv:1311.0851 .
Fan, J., Fan, Y. & Lv, J. (2008), ‘High dimensional covariance matrix estimation using a factor
model’, Journal of Econometrics 147(1), 186–197.
Fan, J., Liao, Y. & Mincheva, M. (2013), ‘Large covariance estimation by thresholding principal
orthogonal complements’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75(4), 603–680.
Friedman, J., Hastie, T. & Tibshirani, R. (2008), ‘Sparse inverse covariance estimation with the
graphical lasso’, Biostatistics 9, 432–441.
Hahsler, M. & Hornik, K. (2017), TSP: Traveling Salesperson Problem (TSP). R package version
1.1-5.
URL: https://CRAN.R-project.org/package=TSP
Hastie, T., Tibshirani, R. & Friedman, J. (2009), The Elements of Statistical Learning; Data
Mining, Inference and Prediction, Second Edition, Springer Verlag, New York.
Huang, J., Liu, N., Pourahmadi, M. & Liu, L. (2006), ‘Covariance matrix selection and estimation
via penalised normal likelihood’, Biometrika 93(1), 85.
Jacob, L., Obozinski, G. & Vert, J.-P. (2009), Group lasso with overlap and graph lasso, in ‘Proceedings of the 26th annual international conference on machine learning’, ACM, pp. 433–440.
Johnstone, I. M. (2001), ‘On the distribution of the largest eigenvalue in principal components
analysis’, Annals of statistics pp. 295–327.
Khare, K., Oh, S.-Y. & Rajaratnam, B. (2014), ‘A convex pseudolikelihood framework for high
dimensional partial correlation estimation with convergence guarantees’, Journal of the Royal
Statistical Society: Series B (Statistical Methodology) .
Lam, C. & Fan, J. (2009), ‘Sparsistency and rates of convergence in large covariance matrix estimation’, Annals of Statistics 37(6B), 4254–4278.
Le Cun, B. B., Denker, J. S., Henderson, D., Howard, R. E., Hubbard, W. & Jackel, L. D. (1990),
Handwritten digit recognition with a back-propagation network, in ‘Advances in neural information processing systems’, Citeseer.
Ledoit, O. & Wolf, M. (2004), ‘A well-conditioned estimator for large-dimensional covariance matrices’, Journal of multivariate analysis 88(2), 365–411.
Levina, E., Rothman, A. & Zhu, J. (2008), ‘Sparse estimation of large covariance matrices via a
nested lasso penalty’, The Annals of Applied Statistics pp. 245–263.
Liu, H., Wang, L. & Zhao, T. (2013), ‘Sparse covariance matrix estimation with eigenvalue constraints’, Journal of Computational and Graphical Statistics .
30

Liu, H., Wang, L. et al. (2017), ‘Tiger: A tuning-insensitive approach for optimally estimating
gaussian graphical models’, Electronic Journal of Statistics 11(1), 241–294.
Luo, X. (2011), ‘Recovering model structures from large low rank and sparse covariance matrix
estimation’, arXiv preprint arXiv:1111.1133 .
Meinshausen, N. & Bühlmann, P. (2006), ‘High dimensional graphs and variable selection with the
lasso’, Annals of Statistics 34, 1436–1462.
Obozinski, G., Jacob, L. & Vert, J. (2011), ‘Group lasso with overlaps: The latent group lasso
approach’, arXiv preprint arXiv:1110.0413 .
Parikh, N. & Boyd, S. P. (2014), ‘Proximal algorithms.’, Foundations and Trends in optimization
1(3), 127–239.
Peng, J., Wang, P., Zhou, N. & Zhu, J. (2009), ‘Partial correlation estimation by joint sparse
regression models’, Journal of the American Statistical Association 104(486).
R Core Team (2017), R: A Language and Environment for Statistical Computing, R Foundation
for Statistical Computing, Vienna, Austria.
URL: https://www.R-project.org/
Rothman, A. (2012), ‘Positive definite estimators of large covariance matrices’, Biometrika
99(3), 733–740.
Rothman, A. J., Levina, E. & Zhu, J. (2009), ‘Generalized thresholding of large covariance matrices’,
Journal of the American Statistical Association 104(485), 177–186.
Rothman, A., Levina, E. & Zhu, J. (2008), ‘Sparse permutation invariant covariance estimation’,
Electronic Journal of Statistics 2, 494–515.
Rothman, A., Levina, E. & Zhu, J. (2010), ‘A new approach to Cholesky-based covariance regularization in high dimensions’, Biometrika 97(3), 539.
Sun, Q., Tan, K. M., Liu, H. & Zhang, T. (2017), ‘Graphical nonconvex optimization for optimal
estimation in gaussian graphical models’, arXiv preprint arXiv:1706.01158 .
Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, J. Royal. Stat. Soc. B.
58, 267–288.
Travers, J. & Milgram, S. (1969), ‘An experimental study of the small world problem’, Sociometry
pp. 425–443.
Tseng, P. (2001), ‘Convergence of a block coordinate descent method for nondifferentiable minimization’, Journal of optimization theory and applications 109(3), 475–494.
Wagaman, A. & Levina, E. (2009), ‘Discovering sparse covariance structures with the isomap’,
Journal of Computational and Graphical Statistics 18(3), 551–572.
Watts, D. J. & Strogatz, S. H. (1998), ‘Collective dynamics of ‘small-world’ networks’, Nature
393(6684), 440–442.
31

Wickham, H. (2009), ggplot2: elegant graphics for data analysis, Springer New York.
URL: http://had.co.nz/ggplot2/book
Won, J.-H., Lim, J., Kim, S.-J. & Rajaratnam, B. (2013), ‘Condition-number-regularized covariance estimation’, Journal of the Royal Statistical Society: Series B (Statistical Methodology)
75(3), 427–450.
Wu, W. B. & Pourahmadi, M. (2003), ‘Nonparametric estimation of large covariance matrices of
longitudinal data’, Biometrika 90(4), 831–844.
Wu, W. B. & Pourahmadi, M. (2009), ‘Banding sample autocovariance matrices of stationary
processes’, Statistica Sinica 19(4), 1755.
Xue, L., Ma, S. & Zou, H. (2012), ‘Positive-definite l1-penalized estimation of large covariance
matrices’, Journal of the American Statistical Association 107(500), 1480–1491.
Yan, X. & Bien, J. (2017), ‘Hierarchical sparse modeling: A choice of two group lasso formulations’,
Statistical Science 32(4), 531–560.
Yu, G. & Bien, J. (2017), ‘Learning local dependence in ordered data’, Journal of Machine Learning
Research 18(42), 1–60.
Yuan, M. (2010), ‘High dimensional inverse covariance matrix estimation via linear programming’,
The Journal of Machine Learning Research 11, 2261–2286.
Yuan, M. & Lin, Y. (2006), ‘Model selection and estimation in regression with grouped variables’,
Journal of the Royal Statistical Society, Series B 68, 49–67.
Yuan, M. & Lin, Y. (2007), ‘Model selection and estimation in the Gaussian graphical model’,
Biometrika 94(1), 19–35.
Zhu, Z. & Liu, Y. (2009), ‘Estimating spatial covariance using penalised likelihood with weighted
l1 penalty’, Journal of Nonparametric Statistics 21(7), 925–942.

32

