High-dimensional covariance estimation by minimizing $\ell_1$-penalized log-determinant divergence

Introduction

The area of high-dimensional statistics deals with estimation in the “large $p$, small $n$” setting, where $p$ and $n$ correspond, respectively, to the dimensionality of the data and the sample size. Such high-dimensional problems arise in a variety of applications, among them remote sensing, computational biology and natural language processing, where the model dimension may be comparable or substantially larger than the sample size. It is well-known that such high-dimensional scaling can lead to dramatic breakdowns in many classical procedures. In the absence of additional model assumptions, it is frequently impossible to obtain consistent procedures when $p\gg n$. Accordingly, an active line of statistical research is based on imposing various restrictions on the model—-for instance, sparsity, manifold structure, or graphical model structure—-and then studying the scaling behavior of different estimators as a function of sample size $n$, ambient dimension $p$ and additional parameters related to these structural assumptions.

Our first result establishes consistency of our estimator $\widehat{\Theta}$ in the elementwise maximum-norm, providing a rate that depends on the tail behavior of the entries in the random matrix $\widehat{\Sigma}^{n}-\Sigma^{*}$. For the special case of sub-Gaussian random vectors with concentration matrices having at most $d$ non-zeros per row, a corollary of our analysis is consistency in spectral norm at rate $|\!|\!|\widehat{\Theta}-\Theta^{*}|\!|\!|_{{2}}={\mathcal{O}}(\sqrt{(d^{2}\,\log p)/n})$, with high probability, thereby strengthening previous results . Under the milder restriction of each element of $X$ having bounded $4m$-th moment, the rate in spectral norm is substantially slower—namely, $|\!|\!|\widehat{\Theta}-\Theta^{*}|\!|\!|_{{2}}={\mathcal{O}}(d\,p^{1/2m}/\sqrt{n})$—highlighting that the familiar logarithmic dependence on the model size $p$ is linked to particular tail behavior of the distribution of $X$. Finally, we show that under the same scalings as above, with probability converging to one, the estimate $\widehat{\Theta}$ correctly specifies the zero pattern of the concentration matrix $\Theta^{*}$.

The remainder of this paper is organized as follows. In Section 2, we set up the problem and give some background. Section 3 is devoted to statements of our main results, as well as discussion of their consequences. Section 4 provides an outline of the proofs, with the more technical details deferred to appendices. In Section 5, we report the results of some simulation studies that illustrate our theoretical predictions.

Background and problem set-up

One motivation for this paper is the problem of Gaussian graphical model selection. A graphical model or a Markov random field is a family of probability distributions for which the conditional independence and factorization properties are captured by a graph. Let $X=(X_{1},X_{2},\ldots,X_{p})$ denote a zero-mean Gaussian random vector; its density can be parameterized by the inverse covariance or concentration matrix $\Theta^{*}=(\Sigma^{*})^{-1}\in\mathcal{S}^{p}_{+}$, and can be written as

corresponds to the problem of Gaussian graphical model selection.

With a slight abuse of notation, we define the sparsity index $s:=|E(\Theta^{*})|$ as the total number of non-zero elements in off-diagonal positions of $\Theta^{*}$; equivalently, this corresponds to twice the number of edges in the case of a Gaussian graphical model. We also define the maximum degree or row cardinality

corresponding to the maximum number of non-zeros in any row of $\Theta^{*}$; this corresponds to the maximum degree in the graph of the underlying Gaussian graphical model. Note that we have included the diagonal entry $\Theta^{*}_{ii}$ in the degree count, corresponding to a self-loop at each vertex.

It is convenient throughout the paper to use graphical terminology, such as degrees and edges, even though the distributional assumptions that we impose, as described in Section 2.3, are milder and hence apply even to distributions that are not Gaussian MRFs.

An important set in this paper is the cone

formed by all symmetric positive semi-definite matrices in $p$ dimensions. We assume that the covariance matrix $\Sigma^{*}$ and concentration matrix $\Theta^{*}$ of the random vector $X$ are strictly positive definite, and so lie in the interior of this cone $\mathcal{S}^{p}_{+}$.

The focus of this paper is a particular type of $M$-estimator for the concentration matrix $\Theta^{*}$, based on minimizing a Bregman divergence between symmetric matrices. A function is of Bregman type if it is strictly convex, continuously differentiable and has bounded level sets . Any such function induces a Bregman divergence of the form $D_{g}(A\|B)=g(A)-g(B)-\left<\nabla g(B),A-B\right>$. From the strict convexity of $g$, it follows that $D_{g}(A\|B)\geq 0$ for all $A$ and $B$, with equality if and only if $A=B$.

As a candidate Bregman function, consider the log-determinant barrier function, defined for any matrix $A\in\mathcal{S}^{p}_{+}$ by

valid for any $A,B\in\mathcal{S}^{p}_{+}$ that are strictly positive definite. This divergence suggests a natural way to estimate concentration matrices—namely, by minimizing the divergence $D_{g}(\Theta^{*}\|\Theta)$—or equivalently, by minimizing the function

In this paper, we analyze a particular instantiation of this strategy. Given $n$ samples, we define the sample covariance matrix

3 Tail conditions

In this section, we describe the tail conditions that underlie our analysis. Since the estimator (11) is based on using the sample covariance $\widehat{\Sigma}^{n}$ as a surrogate for the (unknown) covariance $\Sigma^{*}$, any type of consistency requires bounds on the difference $\widehat{\Sigma}^{n}-\Sigma^{*}$. In particular, we define the following tail condition:

We adopt the convention $1/0:=+\infty$, so that the value $v_{*}=0$ indicates the inequality holds for any $\delta\in(0,\infty)$.

Two important examples of the tail function $f$ are the following:

an exponential-type tail function, meaning that $f(n,\delta)=\exp(c\,n\,\delta^{a})$, for some scalar $c>0$, and exponent $a>0$; and

As might be expected, if $X$ is multivariate Gaussian, then the deviations of sample covariance matrix have an exponential-type tail function with $a=2$. A bit more generally, in the following subsections, we provide broader classes of distributions whose sample covariance entries satisfy exponential and a polynomial tail bounds (see Lemmata 1 and 2 respectively).

Given a larger number of samples $n$, we expect the tail probability bound $1/f(n,\delta)$ to be smaller, or equivalently, for the tail function $f(n,\delta)$ to larger. Accordingly, we require that $f$ is monotonically increasing in $n$, so that for each fixed $\delta>0$, we can define the inverse function

Similarly, we expect that $f$ is monotonically increasing in $\delta$, so that for each fixed $n$, we can define the inverse in the second argument

For future reference, we note a simple consequence of the monotonicity of the tail function $f$—namely

The inverse functions ${\makebox[0.0pt][l]{\hskip 1.29167pt\rule[5.59721pt]{4.47513pt}{0.43057pt}}{n}_{f}}$ and $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}$ play an important role in describing the behavior of our estimator. We provide concrete examples in the following two subsections.

In this subsection, we study the case of i.i.d. observations of sub-Gaussian random variables.

A zero-mean random variable $Z$ is sub-Gaussian if there exists a constant $\sigma\in(0,\infty)$ such that

By the Chernoff bound, this upper bound (16) on the moment-generating function implies a two-sided tail bound of the form

Naturally, any zero-mean Gaussian variable with variance $\sigma^{2}$ satisfies the bounds (16) and (17). In addition to the Gaussian case, the class of sub-Gaussian variates includes any bounded random variable (e.g., Bernoulli, multinomial, uniform), any random variable with strictly log-concave density , and any finite mixture of sub-Gaussian variables.

The following lemma, proved in Appendix D, shows that the entries of the sample covariance based on i.i.d. samples of sub-Gaussian random vector satisfy an exponential-type tail bound with exponent $a=2$. The argument is along the lines of a result due to Bickel and Levina , but with more explicit control of the constants in the error exponent:

Consider a zero-mean random vector $(X_{1},\ldots,X_{p})$ with covariance $\Sigma^{*}$ such that each $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ is sub-Gaussian with parameter $\sigma$. Given $n$ i.i.d. samples, the associated sample covariance $\widehat{\Sigma}^{n}$ satisfies the tail bound

for all \delta\in\big{(}0,\max_{i}(\Sigma^{*}_{ii})\,8(1+4\sigma^{2})\big{)}.

Thus, the sample covariance entries the tail condition $\mathcal{T}(f,v_{*})$ with v_{*}=\big{[}\max_{i}(\Sigma^{*}_{ii})\,8(1+4\sigma^{2})\big{]}^{-1}, and an exponential-type tail function with $a=2$—namely

A little calculation shows that the associated inverse functions take the form

3.2 Tail bounds with moment bounds

In the following lemma, proved in Appendix E, we show that given i.i.d. observations from random variables with bounded moments, the sample covariance entries satisfy a polynomial-type tail bound. See the papers for related results on tail bounds for variables with bounded moments.

For i.i.d. samples $\{X^{(k)}_{i}\}_{k=1}^{n}$, the sample covariance matrix $\widehat{\Sigma}^{n}$ satisfies the bound

Thus, in this case, the sample covariance satisfies the tail condition $\mathcal{T}(f,v_{*})$ with $v_{*}=0$, so that the bound holds for all $\delta\in(0,\infty)$, and with the polynomial-type tail function

Finally, a little calculation shows that in this case, the inverse tail functions take the form

Main results and some consequences

Our results involve some quantities involving the Hessian of the log-determinant barrier (6), evaluated at the true concentration matrix $\Theta^{*}$. Using standard results on matrix derivatives , it can be shown that this Hessian takes the form

We define the set of non-zero off-diagonal entries in the model concentration matrix $\Theta^{*}$:

Our analysis keeps explicit track of these quantities, so that they can scale in a non-trivial manner with the problem dimension $p$.

We assume the Hessian satisfies the following type of mutual incoherence or irrepresentable condition:

There exists some $\alpha\in(0,1]$ such that

The underlying intuition is that this assumption imposes control on the influence that the non-edge terms, indexed by $S^{c}$, can have on the edge-based terms, indexed by $S$. It is worth noting that a similar condition for the Lasso, with the covariance matrix $\Sigma^{*}$ taking the place of the matrix $\Gamma^{*}$ above, is necessary and sufficient for support recovery using the ordinary Lasso . See Section 3.4 for illustration of the form taken by Assumption 1 for specific graphical models.

A remark on notation: although our analysis allows the quantities $K_{\Sigma^{*}},K_{\Gamma^{*}}$ as well as the model size $p$ and maximum node-degree $d$ to grow with the sample size $n$, we suppress this dependence on $n$ in their notation.

In the theorem statement, the choice of regularization constant $\lambda_{n}$ is specified in terms of a user-defined parameter $\tau>2$. Larger choices of $\tau$ yield faster rates of convergence in the probability with which the claims hold, but also lead to more stringent requirements on the sample size.

Consider a distribution satisfying the incoherence assumption (28) with parameter $\alpha\in(0,1]$, and the tail condition (12) with parameters $\mathcal{T}(f,v_{*})$. Let $\widehat{\Theta}$ be the unique optimum of the log-determinant program (11) with regularization parameter $\lambda_{n}=(8/\alpha)\,\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$ for some $\tau>2$. Then, if the sample size is lower bounded as

then with probability greater than $1-1/p^{\tau-2}\rightarrow 1$, we have:

It specifies an edge set $E(\widehat{\Theta})$ that is a subset of the true edge set $E(\Theta^{*})$, and includes all edges $(i,j)$ with |\Theta^{*}_{ij}|>\big{\{}2\big{(}1+8\alpha^{-1}\big{)}K_{\Gamma^{*}}\big{\}}\;\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau}).

We now discuss the consequences of Theorem 1 for distributions in which the sample covariance satisfies an exponential-type tail bound with exponent $a=2$. In particular, recall from Lemma 1 that such a tail bound holds when the variables are sub-Gaussian.

Under the same conditions as Theorem 1, suppose moreover that the variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ are sub-Gaussian with parameter $\sigma$, and the samples are drawn independently. Then if the sample size $n$ satisfies the bound

where C_{1}:=\big{\{}48\sqrt{2}\,(1+4\sigma^{2})\,\max_{i}(\Sigma^{*}_{ii})\,\max\{K_{\Sigma^{*}}K_{\Gamma^{*}},K_{\Sigma^{*}}^{3}K_{\Gamma^{*}}^{2}\}\big{\}}^{2}, then with probability greater than $1-1/p^{\tau-2}$, the estimate $\widehat{\Theta}$ satisfies the bound,

From Lemma 1, when the rescaled variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ are sub-Gaussian with parameter $\sigma$, the sample covariance entries satisfies a tail bound $\mathcal{T}(f,v_{*})$ with with v_{*}=\big{[}\max_{i}(\Sigma^{*}_{ii})\,8(1+4\sigma^{2})\big{]}^{-1} and $f(n,\delta)=(1/4)\exp(c_{*}n\delta^{2})$, where c_{*}=\big{[}128(1+4\sigma^{2})^{2}\max_{i}(\Sigma^{*}_{ii})^{2}\big{]}^{-1}. As a consequence, for this particular model, the inverse functions $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$ and ${\makebox[0.0pt][l]{\hskip 1.29167pt\rule[5.59721pt]{4.47513pt}{0.43057pt}}{n}_{f}}(\delta,p^{\tau})$ take the form

Substituting these forms into the claim of Theorem 1 and doing some simple algebra yields the stated corollary. ∎

2.2 Polynomial-type tails

We now state a corollary for the case of a polynomial-type tail function, such as those ensured by the case of random variables with appropriately bounded moments.

Under the assumptions of Theorem 1, suppose the rescaled variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ have $4m^{th}$ moments upper bounded by $K_{m}$, and the sampling is i.i.d. Then if the sample size $n$ satisfies the bound

where C_{2}:=\big{\{}12m\,[m(K_{m}+1)]^{\frac{1}{2m}}\,\max_{i}(\Sigma^{*}_{ii})\max\{K_{\Sigma^{*}}^{2}K_{\Gamma^{*}},K_{\Sigma^{*}}^{4}K_{\Gamma^{*}}^{2}\}\big{\}}^{2}, then with probability greater than $1-1/p^{\tau-2}$, the estimate $\widehat{\Theta}$ satisfies the bound,

Recall from Lemma 2 that when the rescaled variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ have bounded $4m^{th}$ moments, then the sample covariance $\widehat{\Sigma}$ satisfies the tail condition $\mathcal{T}(f,v_{*})$ with $v_{*}=0$, and with $f(n,\delta)=c_{*}n^{m}\delta^{2m}$ with $c_{*}$ defined as c_{*}=1/\big{\{}m^{2m+1}2^{2m}(\max_{i}\Sigma^{*}_{ii})^{2m}\,(K_{m}+1)\big{\}}. As a consequence, for this particular model, the inverse functions take the form

The claim then follows by substituting these expressions into Theorem 1 and performing some algebra. ∎

3 Model selection consistency

Part (b) of Theorem 1 asserts that the edge set $E(\widehat{\Theta})$ returned by the estimator is contained within the true edge set $E(\Theta^{*})$—meaning that it correctly excludes all non-edges—and that it includes all edges that are “large”, relative to the $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$ decay of the error. The following result, essentially a minor refinement of Theorem 1, provides sufficient conditions linking the sample size $n$ and the minimum value

for model selection consistency. More precisely, define the event

that the estimator $\widehat{\Theta}$ has the same edge set as $\Theta^{*}$, and moreover recovers the correct signs on these edges. With this notation, we have:

Under the same conditions as Theorem 1, suppose that the sample size satisfies the lower bound

Then the estimator is model selection consistent with high probability as $p\rightarrow\infty$,

In comparison to Theorem 1, the sample size requirement (37) differs only in the additional term $\frac{2K_{\Gamma^{*}}(1+\frac{8}{\alpha})}{\theta_{\operatorname{min}}}$ involving the minimum value. This term can be viewed as constraining how quickly the minimum can decay as a function of $(n,p)$, as we illustrate with some concrete tail functions.

Recall the setting of Section 2.3.1, where the random variables $\{X^{(k)}_{i}/\sqrt{\Sigma^{*}_{ii}}\}$ are sub-Gaussian with parameter $\sigma$. Let us suppose that the parameters $(K_{\Gamma^{*}},K_{\Sigma^{*}},\alpha)$ are viewed as constants (not scaling with $(p,d)$. Then, using the expression (32) for the inverse function ${\makebox[0.0pt][l]{\hskip 1.29167pt\rule[5.59721pt]{4.47513pt}{0.43057pt}}{n}_{f}}$ in this setting, a corollary of Theorem 2 is that a sample size

is sufficient for model selection consistency with probability greater than $1-1/p^{\tau-2}$. Alternatively, we can state that $n=\Omega(\tau d^{2}\log p)$ samples are sufficient, as along as the minimum value scales as $\theta_{\operatorname{min}}=\Omega(\sqrt{\frac{\log p}{n}})$.

3.2 Polynomial-type tails

Recall the setting of Section 2.3.2, where the rescaled random variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ have bounded $4m^{th}$ moments. Using the expression (34) for the inverse function ${\makebox[0.0pt][l]{\hskip 1.29167pt\rule[5.59721pt]{4.47513pt}{0.43057pt}}{n}_{f}}$ in this setting, a corollary of Theorem 2 is that a sample size

is sufficient for model selection consistency with probability greater than $1-1/p^{\tau-2}$. Alternatively, we can state than $n=\Omega(d^{2}p^{\tau/m})$ samples are sufficient, as long as the minimum value scales as $\theta_{\operatorname{min}}=\Omega(p^{\tau/(2m)}/{\sqrt{n}})$.

4 Comparison to neighbor-based graphical model selection

For comparison, consider the application of Theorem 2 to the case where the variables are sub-Gaussian (which includes the Gaussian case). For this setting, we have seen that the scaling required by Theorem 2 is $n=\Omega(\{d^{2}+\theta_{\operatorname{min}}^{-2}\}\log p)$, so that the dependence of the log-determinant approach in $\theta_{\operatorname{min}}$ is identical, but it depends quadratically on the maximum degree $d$. We suspect that that the quadratic dependence $d^{2}$ might be an artifact of our analysis, but have not yet been able to reduce it to $d$. Otherwise, the primary difference between the two methods is in the nature of the irrepresentability assumptions that are imposed: our method requires Assumption 1 on the Hessian $\Gamma^{*}$, whereas the neighborhood-based method imposes this same type of condition on a set of $p$ covariance matrices, each of size $(p-1)\times(p-1)$, one for each node of the graph. Below we show two cases where the Lasso irrepresentability condition holds, while the log-determinant requirement fails. However, in general, we do not know whether the log-determinant irrepresentability strictly dominates its analog for the Lasso.

Consider the following Gaussian graphical model example from Meinshausen . Figure 2(a) shows a diamond-shaped graph $G=(V,E)$, with vertex set $V=\{1,2,3,4\}$ and edge-set as the fully connected graph over $V$ with the edge $(1,4)$ removed.

an inequality which holds for all $\rho\in(-0.2017,0.2017)$. Note that the upper value $0.2017$ is just below the necessary threshold discussed by Meinshausen . On the other hand, the irrepresentability condition for the Lasso requires only that $2|\rho|<1$, i.e., $\rho\in(-0.5,0.5)$. Thus, in the regime $|\rho|\in[0.2017,0.5)$, the Lasso irrepresentability condition holds while the log-determinant counterpart fails.

4.2 Illustration of irrepresentability: Star graphs

A second interesting example is the star-shaped graphical model, illustrated in Figure 2(b), which consists of a single hub node connected to the rest of the spoke nodes. We consider a four node graph, with vertex set $V=\{1,2,3,4\}$ and edge-set $E=\{(1,s)\mid s\in\{2,3,4\}\}$. The covariance matrix $\Sigma^{*}$ is parameterized the correlation parameter $\rho\in$: the diagonal entries are set to $\Sigma^{*}_{ii}=1$, for all $i\in V$; the entries corresponding to edges are set to $\Sigma^{*}_{ij}=\rho$ for $(i,j)\in E$; while the non-edge entries are set as $\Sigma^{*}_{ij}=\rho^{2}$ for $(i,j)\notin E$. Consequently, for this particular example, Assumption 1 reduces to the constraint $|\rho|(|\rho|+2)<1$, which holds for all $\rho\in(-0.414,0.414)$. The irrepresentability condition for the Lasso on the other hand allows the full range $\rho\in(-1,1)$. Thus there is again a regime, $|\rho|\in[0.414,1)$, where the Lasso irrepresentability condition holds while the log-determinant counterpart fails.

5 Rates in Frobenius and spectral norm

We now derive some corollaries of Theorem 1 concerning estimation of $\Theta^{*}$ in Frobenius norm, as well as the spectral norm. Recall that $s=|E(\Theta^{*})|$ denotes the total number of off-diagonal non-zeros in $\Theta^{*}$.

Under the same assumptions as Theorem 1, with probability at least $1-1/p^{\tau-2}$, the estimator $\widehat{\Theta}$ satisfies

With the shorthand notation $\nu:=2K_{\Gamma^{*}}(1+8/\alpha)\;\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$, Theorem 1 guarantees that, with probability at least $1-1/p^{\tau-2}$, $\|\widehat{\Theta}-\Theta^{*}\|_{\infty}\leq\nu$. Since the edge set of $\widehat{\Theta}$ is a subset of that of $\Theta^{*}$, and $\Theta^{*}$ has at most $p+s$ non-zeros (including the diagonal), we conclude that

from which the bound (41a) follows. On the other hand, for a symmetric matrix, we have

using the definition of the $\nu_{\infty}$-operator norm, and the fact that $\widehat{\Theta}$ and $\Theta^{*}$ have at most $d$ non-zeros per row. Since the Frobenius norm upper bounds the spectral norm, the bound (41b) follows.

For the exponential tail function case where the rescaled random variables $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ are sub-Gaussian with parameter $\sigma$, we can use the expression (32) for the inverse function $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}$ to derive rates in Frobenius and spectral norms. When the quantities $K_{\Gamma^{*}},K_{\Sigma^{*}},\alpha$ remain constant, these bounds can be summarized succinctly as a sample size $n=\Omega(d^{2}\log p)$ is sufficient to guarantee the bounds

with probability at least $1-1/p^{\tau-2}$.

5.2 Polynomial-type tails

Similarly, let us again consider the polynomial tail case, in which the rescaled variates $X_{i}/\sqrt{\Sigma^{*}_{ii}}$ have bounded $4m^{th}$ moments and the samples are drawn i.i.d. Using the expression (34) for the inverse function we can derive rates in the Frobenius and spectral norms. When the quantities $K_{\Gamma^{*}},K_{\Sigma^{*}},\alpha$ are viewed as constant, we are guaranteed that a sample size $n=\Omega(d^{2}\,p^{\tau/m})$ is sufficient to guarantee the bounds

with probability at least $1-1/p^{\tau-2}$.

6 Rates for the covariance matrix estimate

Finally, we describe some bounds on the estimation of the covariance matrix $\Sigma^{*}$. By Lemma 3, the estimated concentration matrix $\widehat{\Theta}$ is positive definite, and hence can be inverted to obtain an estimate of the covariance matrix, which we denote as $\widehat{\widehat{\Sigma}}:=(\widehat{\Theta})^{-1}$.

Under the same assumptions as Theorem 1, with probability at least $1-1/p^{\tau-2}$, the following bounds hold.

where C_{3}=2K_{\Sigma^{*}}^{2}K_{\Gamma^{*}}\Big{(}1+\frac{8}{\alpha}\Big{)} and C_{4}=6K_{\Sigma^{*}}^{3}K_{\Gamma^{*}}^{2}\Big{(}1+\frac{8}{\alpha}\Big{)}^{2}.

The proof involves certain lemmata and derivations that are parts of the proofs of Theorems 1 and 2, so that we defer it to Section 4.5.

Proofs of main result

In this section, we work through the proofs of Theorems 1 and 2. We break down the proofs into a sequence of lemmas, with some of the more technical aspects deferred to appendices.

Our proofs are based on a technique that we call a primal-dual witness method, used previously in analysis of the Lasso . It involves following a specific sequence of steps to construct a pair $(\widetilde{\Theta},\widetilde{Z})$ of symmetric matrices that together satisfy the optimality conditions associated with the convex program (11) with high probability. Thus, when the constructive procedure succeeds, $\widetilde{\Theta}$ is equal to the unique solution $\widehat{\Theta}$ of the convex program (11), and $\widetilde{Z}$ is an optimal solution to its dual. In this way, the estimator $\widehat{\Theta}$ inherits from $\widetilde{\Theta}$ various optimality properties in terms of its distance to the truth $\Theta^{*}$, and its recovery of the signed sparsity pattern. To be clear, our procedure for constructing $\widetilde{\Theta}$ is not a practical algorithm for solving the log-determinant problem (11), but rather is used as a proof technique for certifying the behavior of the $M$-estimator (11).

The following result is proved in Appendix A:

where $\widehat{Z}$ is an element of the subdifferential $\partial\|\widehat{\Theta}\|_{1,\operatorname{off}}$.

Based on this lemma, we construct the primal-dual witness solution $(\widetilde{\Theta},\widetilde{Z})$ as follows:

We determine the matrix $\widetilde{\Theta}$ by solving the restricted log-determinant problem

Note that by construction, we have $\widetilde{\Theta}\succ 0$, and moreover $\widetilde{\Theta}_{S^{c}}=0$.

We choose $\widetilde{Z}_{S}$ as a member of the sub-differential of the regularizer $\|\cdot\|_{1,\operatorname{off}}$, evaluated at $\widetilde{\Theta}$.

which ensures that constructed matrices $(\widetilde{\Theta},\widetilde{Z})$ satisfy the optimality condition (48).

We verify the strict dual feasibility condition

To clarify the nature of the construction, steps (a) through (c) suffice to obtain a pair $(\widetilde{\Theta},\widetilde{Z})$ that satisfy the optimality conditions (48), but do not guarantee that $\widetilde{Z}$ is an element of sub-differential $\partial\|\widetilde{\Theta}\|_{1,\operatorname{off}}$. By construction, specifically step (b) of the construction ensures that the entries $\widetilde{Z}$ in $S$ satisfy the sub-differential conditions, since $\widetilde{Z}_{S}$ is a member of the sub-differential of $\partial\|\widetilde{\Theta}_{S}\|_{1,\operatorname{off}}$. The purpose of step (d), then, is to verify that the remaining elements of $\widetilde{Z}$ satisfy the necessary conditions to belong to the sub-differential.

In the analysis to follow, some additional notation is useful. We let $W$ denote the “effective noise” in the sample covariance matrix $\widehat{\Sigma}$, namely

Second, we use $\Delta=\widetilde{\Theta}-\Theta^{*}$ to measure the discrepancy between the primal witness matrix $\widetilde{\Theta}$ and the truth $\Theta^{*}$. Finally, recall the log-determinant barrier $g$ from equation (6). We let $R(\Delta)$ denote the difference of the gradient $\nabla g(\widetilde{\Theta})={\widetilde{\Theta}}^{-1}$ from its first-order Taylor expansion around $\Theta^{*}$. Using known results on the first and second derivatives of the log-determinant function (see p. 641 in Boyd and Vandenberghe ), this remainder takes the form

2 Auxiliary results

We begin by stating and proving a lemma that provides sufficient (deterministic) conditions for strict dual feasibility to hold, so that $\|\widetilde{Z}_{S^{c}}\|_{\infty}<1$.

Then the matrix $\widetilde{Z}_{S^{c}}$ constructed in step (c) satisfies $\|\widetilde{Z}_{S^{c}}\|_{\infty}<1$, and therefore $\widetilde{\Theta}=\widehat{\Theta}$.

Using the definitions (51) and (52), we can re-write the stationary condition (48) in an alternative but equivalent form

In terms of the disjoint decomposition $S$ and $S^{c}$, equation (54) can be re-written as two blocks of linear equations as follows:

Here we have used the fact that $\Delta_{S^{c}}=0$ by construction.

Since $\Gamma^{*}_{SS}$ is invertible, we can solve for $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.96916pt}{0.43057pt}}{\Delta}_{S}$ from equation (55a) as follows:

Substituting this expression into equation (55b), we can solve for $\widetilde{Z}_{S^{c}}$ as follows:

Recalling Assumption 1—namely, that |\!|\!|\Gamma^{*}_{S^{c}S}\big{(}\Gamma^{*}_{SS}\big{)}^{-1}|\!|\!|_{{\infty}}\leq(1-\alpha)—we have

where we have used the fact that $\|\makebox[0.0pt][l]{\hskip 2.16669pt\rule[8.5139pt]{3.21942pt}{0.43057pt}}{\widetilde{Z}}_{S}\|_{\infty}\leq 1$, since $\widetilde{Z}$ belongs to the sub-differential of the norm $\|\cdot\|_{1,\operatorname{off}}$ by construction. Finally, applying assumption (53) from the lemma statement, we have

2.2 Control of remainder term

Our next step is to relate the behavior of the remainder term (52) to the deviation $\Delta=\widetilde{\Theta}-\Theta^{*}$.

We provide the proof of this lemma in Appendix B. It is straightforward, based on standard matrix expansion techniques.

2.4 Sufficient conditions for sign consistency

We now show how a lower bound on the minimum value $\theta_{\operatorname{min}}$, when combined with Lemma 6, allows us to guarantee sign consistency of the primal witness matrix $\widetilde{\Theta}_{S}$.

Suppose the minimum absolute value $\theta_{\operatorname{min}}$ of non-zero entries in the true concentration matrix $\Theta^{*}$ is lower bounded as

then $\textrm{sign}(\widetilde{\Theta}_{S})=\textrm{sign}(\Theta^{*}_{S})$ holds.

This claim follows from the bound (62) combined with the bound (60) ,which together imply that for all $(i,j)\in S$, the estimate $\widetilde{\Theta}_{ij}$ cannot differ enough from $\Theta^{*}_{ij}$ to change sign.

2.5 Control of noise term

The final ingredient required for the proofs of Theorems 1 and 2 is control on the sampling noise $W=\widehat{\Sigma}-\Sigma^{*}$. This control is specified in terms of the decay function $f$ from equation (12).

For any $\tau>2$ and sample size $n$ such that $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})\leq 1/v_{*}$, we have

Using the definition (12) of the decay function $f$, and applying the union bound over all $p^{2}$ entries of the noise matrix, we obtain that for all $\delta\leq 1/v_{*}$,

Setting $\delta=\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$ yields that

as claimed. Here the last equality follows since $f(n,\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau}))=p^{\tau}$, using the definition (14) of the inverse function $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}$. ∎

3 Proof of Theorem 1

We proceed by verifying that assumption (53) of Lemma 4 holds. Recalling the choice of regularization penalty $\lambda_{n}=(8/\alpha)\,\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$, we have $\|W\|_{\infty}\leq(\alpha/8)\lambda_{n}$. In order to establish condition (53) it remains to establish the bound $\|R(\Delta)\|_{\infty}\leq\frac{\alpha\,\lambda_{n}}{8}$. We do so in two steps, by using Lemmas 6 and 5 consecutively. First, we show that the precondition (59) required for Lemma 6 to hold is satisfied under the specified conditions on $n$ and $\lambda_{n}$. From Lemma 8 and our choice of regularization constant $\lambda_{n}=(8/\alpha)\,\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})$,

provided $\makebox[0.0pt][l]{\hskip 2.08334pt\rule[8.23611pt]{2.24303pt}{0.43057pt}}{\delta}_{f}(n,p^{\tau})\leq 1/v_{*}$. From the lower bound (29) and the monotonicity (15) of the tail inverse functions, we have

showing that the assumptions of Lemma 6 are satisfied. Applying this lemma, we conclude that

Turning next to Lemma 5, we see that its assumption $\|\Delta\|_{\infty}\leq\frac{1}{3\,K_{\Sigma^{*}}d}$ holds, by applying equations (64) and (65). Consequently, we have

as required, where the final inequality follows from our condition (29) on the sample size, and the monotonicity property (15).

4 Proof of Theorem 2

We now turn to the proof of Theorem 2. A little calculation shows that the assumed lower bound (37) on the sample size $n$ and the monotonicity property (15) together guarantee that

Proceeding as in the proof of Theorem 1, with probability at least $1-1/p^{\tau-2}$, we have the equality $\widetilde{\Theta}=\widehat{\Theta}$, and also that $\|\widetilde{\Theta}-\Theta^{*}\|_{\infty}\leq\theta_{\operatorname{min}}/2$. Consequently, Lemma 7 can be applied, guaranteeing that $\textrm{sign}(\Theta^{*}_{ij})=\textrm{sign}(\widetilde{\Theta}_{ij})$ for all $(i,j)\in E$. Overall, we conclude that with probability at least $1-1/p^{\tau-2}$, the sign consistency condition $\textrm{sign}(\Theta^{*}_{ij})=\textrm{sign}(\widehat{\Theta}_{ij})$ holds for all $(i,j)\in E$, as claimed.

5 Proof of Corollary 4

With the shorthand $\widehat{\Delta}=\widehat{\Theta}-\Theta^{*}$, we have

From the definition (52) of the residual $R(\cdot)$, this difference can be written as

where J:=\sum_{k=0}^{\infty}(-1)^{k}\big{(}{{\Theta}^{*}}^{-1}\widehat{\Delta}\big{)}^{k} has norm $|\!|\!|J^{T}|\!|\!|_{{\infty}}\leq 3/2$.

The quantity $L(\widehat{\Delta})$ in turn can be bounded as follows,

where we used the inequality that $\|\widehat{\Delta}u\|_{\infty}\leq\|\widehat{\Delta}\|_{\infty}\|u\|_{1}$. Simplifying further, we obtain

where we have used the fact that $|\!|\!|{\Theta^{*}}^{-1}|\!|\!|_{{1}}=|\!|\!|[{\Theta^{*}}^{-1}]^{T}|\!|\!|_{{\infty}}=|\!|\!|{\Theta^{*}}^{-1}|\!|\!|_{{\infty}}$, which follows from the symmetry of ${\Theta^{*}}^{-1}$. Combining the pieces, we obtain

where the last inequality uses the bound $|\!|\!|J|\!|\!|_{{\infty}}\leq 3/2$. (Proceeding as in the proof of Lemma 5, this bound holds conditioned on $\mathcal{A}$, and for the sample size specified in the theorem statement.) In turn, the term $L(\widehat{\Delta})$ can be bounded as

Experiments

Figure 3 illustrates the three types of graphs used in our simulations: chain graphs (panel (a)), four-nearest neighbor lattices or grids (panel (b)), and star-shaped graphs (panel (c)). For the chain and grid graphs, the maximal node degree $d$ is fixed by definition, to $d=2$ for chains, and $d=4$ for the grids. Consequently, these graphs can capture the dependence of the required sample size $n$ only as a function of the graph size $p$, and the parameters $(K_{\Sigma^{*}}$, $K_{\Gamma^{*}}$, $\theta_{\operatorname{min}}$). The star graph allows us to vary both $d$ and $p$, since the degree of the central hub can be varied between $1$ and $p-1$. For each graph type, we varied the size of the graph $p$ in different ranges, from $p=64$ upwards to $p=375$.

For the chain and star graphs, we define a covariance matrix $\Sigma^{*}$ with entries $\Sigma^{*}_{ii}=1$ for all $i=1,\ldots,p$, and $\Sigma^{*}_{ij}=\rho$ for all $(i,j)\in E$ for specific values of $\rho$ specified below. Note that these covariance matrices are sufficient to specify the full model. For the four-nearest neighbor grid graph, we set the entries of the concentration matrix $\Theta^{*}_{ij}=\omega$ for $(i,j)\in E$, with the value $\omega$ specified below. In all cases, we set the regularization parameter $\lambda_{n}$ proportional to $\sqrt{\log(p)/n}$, as suggested by Theorems 1 and 2, which is reasonable since the main purpose of these simulations is to illustrate our theoretical results. However, for general data sets, the relevant theoretical parameters cannot be computed (since the true model is unknown), so that a data-driven approach such as cross-validation might be required for selecting the regularization parameter $\lambda_{n}$.

Given a Gaussian graphical model instance, and the number of samples $n$, we drew $N=100$ batches of $n$ independent samples from the associated multivariate Gaussian distribution. We estimated the probability of correct model selection as the fraction of the $N=100$ trials in which the estimator recovers the signed-edge set exactly.

Then, as a corollary of Theorem 2, a sample size of order

is sufficient for model selection consistency with probability greater than $1-1/p^{\tau-2}$. In the subsections to follow, we investigate how the empirical sample size $n$ required for model selection consistency scales in terms of graph size $p$, maximum degree $d$, as well as the “model-complexity” term $K$ defined above.

Panel (a) of Figure 4 plots the probability of correct signed edge-set recovery against the sample size $n$ for a chain-structured graph of three different sizes. For these chain graphs, regardless of the number of nodes $p$, the maximum node degree is constant $d=2$, while the edge covariances are set as $\Sigma_{ij}=0.2$ for all $(i,j)\in E$, so that the quantities $(K_{\Sigma^{*}},K_{\Gamma^{*}},\alpha)$ remain constant. Each of the curve in panel (a) corresponds to a different graph size $p$. For each curve, the probability of success starts at zero (for small sample sizes $n$), but then transitions to one as the sample size is increased. As would be expected, it is more difficult to perform model selection for larger graph sizes, so that (for instance) the curve for $p=375$ is shifted to the right relative to the curve for $p=64$. Panel (b) of Figure 4 replots the same data, with the horizontal axis rescaled by $(1/\log p)$. This scaling was chosen because for sub-Gaussian tails, our theory predicts that the sample size should scale logarithmically with $p$ (see equation (71)). Consistent with this prediction, when plotted against the rescaled sample size $n/\log p$, the curves in panel (b) all stack up. Consequently, the ratio $(n/\log p)$ acts as an effective sample size in controlling the success of model selection, consistent with the predictions of Theorem 2 for sub-Gaussian variables.

Figure 5 shows the same types of plots for a star-shaped graph with fixed maximum node degree $d=40$, and Figure 6 shows the analogous plots for a grid graph with fixed degree $d=4$. As in the chain case, these plots show the same type of stacking effect in terms of the scaled sample size $n/\log p$, when the degree $d$ and other parameters ($(\alpha,K_{\Gamma^{*}},K_{\Sigma^{*}})$) are held fixed.

2 Dependence on the maximum node degree

Panel (a) of Figure 7 plots the probability of correct signed edge-set recovery against the sample size $n$ for star-shaped graphs; each curve corresponds to a different choice of maximum node degree $d$, allowing us to investigate the dependence of the sample size on this parameter. So as to control these comparisons, the models are chosen such that quantities other than the maximum node-degree $d$ are fixed: in particular, we fix the number of nodes $p=200$, and the edge covariance entries are set as $\Sigma^{*}_{ij}=2.5/d$ for $(i,j)\in E$ so that the quantities $(K_{\Sigma^{*}},K_{\Gamma^{*}},\alpha)$ remain constant. The minimum value $\theta_{\operatorname{min}}$ in turn scales as $1/d$. Observe how the plots in panel (a) shift to the right as the maximum node degree $d$ is increased, showing that star-shaped graphs with higher degrees are more difficult. In panel (b) of Figure 7, we plot the same data versus the rescaled sample size $n/d$. Recall that if all the curves were to stack up under this rescaling, then it means the required sample size $n$ scales linearly with $d$. These plots are closer to aligning than the unrescaled plots, but the agreement is not perfect. In particular, observe that the curve $d$ (right-most in panel (a)) remains a bit to the right in panel (b), which suggests that a somewhat more aggressive rescaling—perhaps $n/d^{\gamma}$ for some $\gamma\in(1,2)$—is appropriate.

Note that for $\theta_{\operatorname{min}}$ scaling as $1/d$, the sufficient condition from Theorem 2, as summarized in equation (71), is $n=\Omega(d^{2}\log p)$, which appears to be overly conservative based on these data. Thus, it might be possible to tighten our theory under certain regimes.

3 Dependence on covariance and Hessian terms

Next, we study the dependence of the sample size required for model selection consistency on the model complexity term $K$ defined in (70), which is a collection of the quantities $K_{\Sigma^{*}}$, $K_{\Gamma^{*}}$ and $\alpha$ defined by the covariance matrix and Hessian, as well as the minimum value $\theta_{\operatorname{min}}$. Figure 8 plots the probability of correct signed edge-set recovery versus the sample size $n$ for chain graphs. Here each curve corresponds to a different setting of the model complexity factor $K$, but with a fixed number of nodes $p=120$, and maximum node-degree $d=2$. We varied the actor $K$ by varying the value $\rho$ of the edge covariances $\Sigma_{ij}=\rho,\,(i,j)\in E$. Notice how the curves, each of which corresponds to a different model complexity factor, shift rightwards as $K$ is increased so that models with larger values of $K$ require greater number of samples $n$ to achieve the same probability of correct model selection. These rightward-shifts are in qualitative agreement with the prediction of Theorem 1, but we suspect that our analysis is not sharp enough to make accurate quantitative predictions regarding this scaling.

Discussion

Our main results relate the i.i.d. sample size $n$ to various parameters of the problem required to achieve consistency. In addition to the dependence on matrix size $p$, number of edges $s$ and graph degree $d$, our analysis also illustrates the role of other quantities, related to the structure of the covariance matrix $\Sigma^{*}$ and the Hessian of the objective function, that have an influence on consistency rates. Our main assumption is an irrepresentability or mutual incoherence condition, similar to that required for model selection consistency of the Lasso, but involving the Hessian of the log-determinant objective function (11), evaluated at the true model. When the distribution of $X$ is multivariate Gaussian, this Hessian is the Fisher information matrix of the model, and thus can be viewed as an edge-based counterpart to the usual node-based covariance matrix We report some examples where irrepresentability condition for the Lasso hold and the log-determinant condition fails, but we do not know in general if one requirement dominates the other. In addition to these theoretical results, we provided a number of simulation studies showing how the sample size required for consistency scales with problem size, node degrees, and the other complexity parameters identified in our analysis.

There are various interesting questions and possible extensions to this paper. First, in the current paper, we have only derived sufficient conditions for model selection consistency. As in past work on the Lasso , it would also be interesting to derive a converse result—namely, to prove that if the sample size $n$ is smaller than some function of $(p,d,s)$ and other complexity parameters, then regardless of the choice of regularization constant, the log-determinant method fails to recover the correct graph structure. Second, while this paper studies the problem of estimating a fixed graph or concentration matrix, a natural extension would allow the graph to vary over time, a problem setting which includes the case where the observations are dependent. For instance, Zhou et al. study the estimation of the covariance matrix of a Gaussian distribution in a time-varying setting, and it would be interesting to extend results of this paper to this more general setting.

We thank Shuheng Zhou for helpful comments on an earlier draft of this work. Work was partially supported by NSF grant DMS-0605165. Yu also acknowledges support from ARO W911NF-05-1-0104, NSFC-60628102, and a grant from MSRA.

Appendix A Proof of Lemma 3

In this appendix, we show that the regularized log-determinant program (11) has a unique solution whenever $\lambda_{n}>0$, and the diagonal of the sample covariance $\widehat{\Sigma}^{n}$ is strictly positive. By the strict convexity of the log-determinant barrier , if the minimum is attained, then it is unique, so that it remains to show that the minimum is achieved. If $\lambda_{n}>0$, then by Lagrangian duality, the problem can be written in an equivalent constrained form:

As long as $\widehat{\Sigma}^{n}_{ii}>0$ for each $i=1,\ldots,p$, this function is coercive, meaning that it diverges to infinity for any sequence $\|(\Theta^{t}_{11},\ldots,\Theta^{t}_{pp})\|\rightarrow+\infty$. Consequently, the minimum is attained.

Returning to the penalized form (11), by standard optimality conditions for convex programs, a matrix $\widehat{\Sigma}\succ 0$ is optimal if and only belongs to the sub-differential of the objective, or equivalently if and only if there exists a matrix $\widehat{Z}$ in the sub-differential of the off-diagonal norm $\|\cdot\|_{1,\operatorname{off}}$ such that

Appendix B Proof of Lemma 5

By sub-multiplicativity of the $|\!|\!|\cdot|\!|\!|_{{\infty}}$ matrix norm, for any two $p\times p$ matrices $A,B$, we have $|\!|\!|A\,B|\!|\!|_{{\infty}}\leq|\!|\!|A|\!|\!|_{{\infty}}|\!|\!|B|\!|\!|_{{\infty}}$, so that

where we have used the definition of $K_{\Sigma^{*}}$, the fact that $\Delta$ has at most $d$ non-zeros per row/column, and our assumption $\|\Delta\|_{\infty}\;<1/(3K_{\Sigma^{*}})$. Consequently, we have the convergent matrix expansion

where J=\sum_{k=0}^{\infty}(-1)^{k}\big{(}{{\Theta}^{*}}^{-1}\Delta\big{)}^{k}.

We now prove the bound (58) on the remainder as follows. Let $e_{i}$ denote the unit vector with $1$ in position $i$ and zeroes elsewhere. From equation (57), we have

Recall that J=\sum_{k=0}^{\infty}(-1)^{k}\big{(}{{\Theta}^{*}}^{-1}\Delta\big{)}^{k}. By sub-multiplicativity of $|\!|\!|\cdot|\!|\!|_{{\infty}}$ matrix norm, we have

since $|\!|\!|{\Theta^{*}}^{-1}|\!|\!|_{{\infty}}|\!|\!|\Delta|\!|\!|_{{\infty}}<1/3$ from equation (73). Substituting this in (B), we obtain

where the final line follows since $|\!|\!|\Delta|\!|\!|_{{\infty}}\leq d\|\Delta\|_{\infty}$, and since $\Delta$ has at most $d$ non-zeroes per row/column.

Appendix C Proof of Lemma 6

By following the same argument as in Appendix A, we conclude that the restricted problem (49) has a unique optimum $\widetilde{\Theta}$. Let $\widetilde{Z}$ be any member of the sub-differential of $\|\cdot\|_{1,\operatorname{off}}$, evaluated at $\widetilde{\Theta}$. By Lagrangian theory, the witness $\widetilde{\Theta}$ must be an optimum of the associated Lagrangian problem

In fact, since this Lagrangian is strictly convex, $\widetilde{\Theta}$ is the only optimum of this problem. Since the log-determinant barrier diverges as $\Theta$ approaches the boundary of the positive semi-definite cone, we must have $\widetilde{\Theta}\succ 0$. If we take partial derivatives of the Lagrangian with respect to the unconstrained elements $\Theta_{S}$, these partial derivatives must vanish at the optimum, meaning that we have the zero-gradient condition

To be clear, $\Theta$ is the $p\times p$ matrix with entries in $S$ equal to $\Theta_{S}$ and entries in $S^{c}$ equal to zero. Since this zero-gradient condition is necessary and sufficient for an optimum of the Lagrangian problem, it has a unique solution (namely, $\widetilde{\Theta}_{S}$).

where $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.52184pt}{0.43057pt}}{G}$ denotes the vectorized form of $G$. Note that by construction, $F(\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.96916pt}{0.43057pt}}{\Delta}_{S})=\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.96916pt}{0.43057pt}}{\Delta}_{S}$ holds if and only if $G(\Theta^{*}_{S}+\Delta_{S})=G(\Theta_{S})=0$.

where we have used the definition $W=\widehat{\Sigma}-\Sigma^{*}$.

By the definition (61) of the radius $r$, and the assumed upper bound (59), we have $\|\Delta\|_{\infty}\leq r\leq\frac{1}{3K_{\Sigma^{*}}d}$, so that the results of Lemma 5 apply. By using the definition (52) of the remainder, taking the vectorized form of the expansion (57), and restricting to entries in $S$, we obtain the expansion

Using this expansion (79) combined with the expression (77) for $G$, we have

The second term is easy to deal with: using the definition $K_{\Gamma^{*}}=|\!|\!|(\Gamma^{*}_{SS})^{-1}|\!|\!|_{{\infty}}$, we have \|T_{2}\|_{\infty}\leq K_{\Gamma^{*}}\big{(}\|W\|_{\infty}+\lambda_{n}\big{)}\;=\;r/2. It now remains to show that $\|T_{1}\|_{\infty}\leq r/2$. We have

where we used the expanded form (57) of the remainder, Applying the bound (58) from Lemma 5, we obtain

Since $r\leq\frac{1}{3K_{\Sigma^{*}}^{3}K_{\Gamma^{*}}d}$ by assumption (59), we conclude that

Appendix D Proof of Lemma 1

For each pair $(i,j)$ and $\nu>0$, define the event

As the sub-Gaussian assumption is imposed on the variables $\{X^{(k)}_{i}\}$ directly, as in Lemma A.3 of Bickel and Levina , our proof proceeds by first decoupling the products $X^{(k)}_{i}X^{(k)}_{j}$. For each pair $(i,j)$, we define $\rho^{*}_{ij}=\Sigma^{*}_{ij}/\sqrt{\Sigma^{*}_{ii}\Sigma^{*}_{jj}}$, and the rescaled random variables $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{6.66843pt}{0.43057pt}}{X}^{(k)}_{i}:=X^{(k)}_{ij}/\sqrt{\Sigma^{*}_{ii}}$. Noting that the strict positive definiteness of $\Sigma^{*}$ implies that $|\rho^{*}_{ij}|<1$, we can also define the auxiliary random variables

Suppose that each $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{6.66843pt}{0.43057pt}}{X}^{(k)}_{i}$ is sub-Gaussian with parameter $\sigma$. Then for each node pair $(i,j)$, the following properties hold:

For all $k=1,\ldots,n$, the random variables $U^{(k)}_{ij}$ and $V^{(k)}_{ij}$ are sub-Gaussian with parameters $2\sigma$.

where we have used the Cauchy-Schwarz inequality. Since the variables $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{6.66843pt}{0.43057pt}}{X}^{(k)}_{i}$ and $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{6.66843pt}{0.43057pt}}{X}^{(k)}_{j}$ are sub-Gaussian with parameter $\sigma$, we have

so that $U^{(k)}_{ij}$ is sub-Gaussian with parameter $2\sigma$ as claimed. (b) By straightforward algebra, we have the decomposition

which completes the proof of Lemma 9(b). ∎

It remains to control the terms $\sum_{k=1}^{n}(U^{(k)}_{ij})^{2}$ and $\sum_{k=1}^{n}(V^{(k)}_{ij})^{2}$. We do so by exploiting tail bounds for sub-exponential random variables. A zero-mean random variable $Z$ is said to be sub-exponential if there exists a constant $\gamma\in(0,\infty)$ and $\phi\in(0,\infty]$ such that

Note that for $\phi<+\infty$, this requirement is a weakening of sub-Gaussianity, since the inequality is only required to hold on the interval $(-\phi,+\phi)$.

Now consider the variates $Z_{k;ij}:=(U^{(k)}_{ij})^{2}-2\,(1+\rho^{*}_{ij})$. Note that they are zero-mean; we also claim they are sub-exponential.

For all $k\in\{1,\ldots,n\}$ and node-pairs $(i,j)\in V\times V$, the variables

are sub-exponential with parameter $\gamma_{U}=16(1+4\sigma^{2})$ in the interval $(-\phi_{U},\phi_{U})$, with $\phi_{U}=1/(16(1+4\sigma^{2}))$.

for all $\nu\leq 8(\max_{i}\Sigma^{*}_{ii})\,(1+4\sigma^{2})$. A similar argument yields the same tail bound for the deviation involving $V^{(k)}_{ij}$. Consequently, using Lemma 9(b), we conclude that

valid for $\nu\leq 8(\max_{i}\Sigma^{*}_{ii})\,(1+4\sigma^{2})$, as required. It only remains to prove Lemma 10.

it then follows (Thm. 3.2, ) that $Z_{k;ij}$ is sub-exponential with parameter $2B$ in the interval $(-\frac{1}{2B},+\frac{1}{2B})$. We obtain such a bound $B$ as follows. Using the inequality $(a+b)^{m}\leq 2^{m}(a^{m}+b^{m})$, valid for any real numbers $a,b$, we have

where we have used the fact that $|\rho^{*}_{ij}|\leq 1$. The claim of the lemma thus follows. ∎

Appendix E Proof of Lemma 2

Define the random variables $W^{(k)}_{ij}=X^{(k)}_{i}X^{(k)}_{j}-\Sigma^{*}_{ij}$, and note that they have mean zero. By applying the Chebyshev inequality, we obtain

Letting $\mathcal{A}=\{(a_{1},\ldots,a_{n})\,\mid\,a_{i}\in\{0,\ldots,2m\},\,\sum_{i=1}^{n}a_{i}=2m\}$, by the multinomial theorem, we have

where the final equality uses linearity of expectation, and the independence of the variables $\{W^{(k)}_{ij}\}_{k=1}^{n}$.

The quantity $T_{1}$ is equal to the number of ways to put $2m$ balls in $n$ bins such that if a bin contains a ball, it should have at least two balls. Note that this implies there can then be at most $m$ bins containing a ball. Consequently, the term $T_{1}$ is bounded above by the product of the number of ways in which we can choose $m$ out of $n$ bins, and the number of ways in which we can put $2m$ balls into $m$ bins—viz.

Using this inequality, for any $a\in\mathcal{A}_{-1}$, we have

Substituting our bounds on $T_{1}$ and $T_{2}$ into equation (E), we obtain

It thus remains to bound the moments of $W^{(k)}_{ij}$. We have

where we have used the inequality $(a+b)^{2m}\;\leq\;2^{2m}(a^{2m}+b^{2m})$, valid for all real numbers $a$ and $b$. An application of the Cauchy-Schwarz inequality yields

Substituting back into equation (84) yields

Noting that $\Sigma^{*}_{ij}$, $\Sigma^{*}_{jj}$, and $\Sigma^{*}_{ij}$ are all bounded above by $\max_{i}\Sigma^{*}_{ii}$, we obtain

References