Asymptotic normality and optimalities in estimation of large Gaussian graphical models

Introduction

The Gaussian graphical model, a powerful tool for investigating the relationship among a large number of random variables in a complex system, is used in a wide range of scientific applications. A central question for Gaussian graphical models is how to recover the structure of an undirected Gaussian graph. Let $G=(V,E)$ be an undirected graph representing the conditional dependence relationship between components of a random vector $Z=(Z_{1},\ldots,Z_{p})^{T}$ as follows. The vertex set $V=\{V_{1},\ldots,V_{p}\}$ represents the components of $Z$. The edge set $E$ consists of pairs $(i,j)$ indicating the conditional dependence between $Z_{i}$ and $Z_{j}$ given all other components. In applications, the following question is fundamental: Is there an edge between $V_{i}$ and $V_{j}$? It is well known that recovering the structure of an undirected Gaussian graph $G=(V,E)$ is equivalent to recovering the support of the population precision matrix of the data in the Gaussian graphical model. Let

where $\Sigma=(\sigma_{ij})$ is the population covariance matrix. The precision matrix, denoted by $\Omega=(\omega_{ij})$, is defined as the inverse of covariance matrix, $\Omega=\Sigma^{-1}$. There is an edge between $V_{i}$ and $V_{j}$, that is, $(i,j)\in E$, if and only if $\omega_{ij}\neq 0$; see, for example, Lauritzen 1996. Consequently, the support recovery of the precision matrix $\Omega$ yields the recovery of the structure of the graph $G$.

Suppose $n$ i.i.d. $p$-variate random vectors $X^{(1)},X^{(2)},\ldots,X^{(n)}$ are observed from the same distribution as $Z$, that is, the Gaussian $\mathcal{N}(\mu,\Omega^{-1})$. Assume without loss of generality that $\mu=0$ hereafter. In this paper, we address the following two fundamental questions: When is it possible to make statistical inference for each individual entry of a precision matrix $\Omega$ at the parametric $n^{-1/2}$ rate? When and in what sense is it possible to recover the support of $\Omega$ in the presence of some small nonzero $|\omega_{ij}|$?

In spite of an extensive literature on the topic, the fundamental limit of support recovery in the Gaussian graphical model is still largely unknown, let alone an adaptive procedure to achieve the limit.

The methodology we are proposing is a novel regression approach briefly described in Sun and Zhang 2012b. In this regression approach, the main task is not to estimate the slope, as seen in Meinshausen and Bühlmann 2006, Yuan 2010, Cai, Liu and Luo 2011, Cai, Liu and Zhou 2012 and Sun and Zhang 2012a, but to estimate the noise level. For a vector $Z$ of length $p$ and any index subset $A$ of $\{1,2,\ldots,p\}$, we denote by $Z_{A}$ the sub-vector of $Z$ with elements indexed by $A$. Similarly for a matrix $U$ and two index subsets $A$ and $B$ of $\{1,2,\ldots,p\}$, we denote by $U_{A,B}$ the $|A|\times|B|$ sub-matrix of $U$ with elements in rows in $A$ and columns in $B$. Consider $A=\{i,j\}$ with $i\neq j$, so that $Z_{A}=(Z_{i},Z_{j})^{T}$ and $\Omega_{A,A}=\bigl({{\omega_{ii}\atop\omega_{ji}}\enskip{\omega_{ij}\atop\omega_{jj}}}\bigr)$. It is well known that

This observation motivates us to consider the estimation of individual entries of $\Omega$, $\omega_{ii}$ and $\omega_{ij}$, by estimating the noise level in the regression of the two response variables in $A$ against the variables in $A^{c}$. The noise level $\Omega_{A,A}^{-1}$ has only three parameters. When $\Omega$ is sufficiently sparse, a penalized regression approach is proposed in Section 2 to obtain an asymptotically efficient estimation of $\omega_{ij}$ in the following sense: The estimator is asymptotically normal, and its asymptotic variance matches that of the maximum likelihood estimator in the classical setting where the dimension $p$ is a fixed constant. Consider the class of parameter spaces modeling sparse precision matrices with at most $k_{n,p}$ nonzero elements in each column,

where $1\{\cdot\}$ is the indicator function, and $M$ is some constant greater than $1$. The following theorem shows that a necessary and sufficient condition to obtain a $n^{-1/2}$-consistent estimation of $\omega_{ij}$ is $k_{n,p}=O(\frac{\sqrt{n}}{\log p})$, and when $k_{n,p}=o(\frac{\sqrt{n}}{\log p})$ the procedure to be proposed in Section 2 is asymptotically efficient.

Let $X^{(i)}{\stackrel{{\scriptstyle\mathit{i.i.d.}}}{{\sim}}}\mathcal{N}_{p}(\mu,\Sigma)$, $i=1,2,\ldots,n$. Assume that $3\leq k_{n,p}\leq c_{0}n/\log p$ with a sufficiently small constant $c_{0}>0$ and $p\geq k_{n,p}^{\nu}$ with some $\nu>2$.

There exists a constant $\varepsilon_{0}>0$ such that

Moreover, the minimax risk of estimating $\omega_{ij}$ over the class $\mathcal{G}_{0}(M,k_{n,p})$ satisfies

uniformly in $(i,j)$, provided that $n=O(p^{\xi})$ with some $\xi>0$.

The estimator $\hat{\omega}_{ij}$ defined in (10) in Section 2 is rate optimal in the sense of

Furthermore, the estimator $\hat{\omega}_{ij}$ is asymptotically efficient when $k_{n,p}=o(\frac{\sqrt{n}}{\log p})$, that is, with $F_{ij}=(\omega_{ii}\omega_{jj}+\omega_{ij}^{2})^{-1}$ being the Fisher information for estimating $\omega_{ij}$ and ${\hat{F}}_{ij}=({\hat{\omega}}_{ii}{\hat{\omega}}_{jj}+{\hat{\omega}}_{ij}^{2})^{-1}$ its estimate,

The lower bound is established through Le Cam’s lemma and a novel construction of a subset of sparse precision matrices. An important implication of the lower bound is that the difficulty of support recovery for sparse precision matrices is different from that for sparse covariance matrices when

$k_{n,p}\gg(\frac{\sqrt{n}}{\log p})$, and when $k_{n,p}=o(\frac{\sqrt{n}}{\log p})$ the difficulty of support recovery for sparse precision matrices is just the same as that for sparse covariance matrices.

The proposed estimator was briefly described in Sun and Zhang 2012b along with a statement of the efficiency of the estimator without proof under the sparsity assumption $k_{n,p}=o(n^{-1/2}\log p)$. While we are working on the delicate issue of the necessity of the sparsity condition $k_{n,p}=o(n^{1/2}/\log p)$ and the optimality of the method for support recovery and estimation under the general sparsity condition $k_{n,p}=o(n/\log p)$, Liu 2013 developed $p$-values for testing $\omega_{ij}=0$ and related FDR control methods under the stronger sparsity condition $k_{n,p}=o(n^{1/2}/\log p)$. However, his method cannot be directly converted into confidence intervals, and the optimality of his method is unclear under either sparsity conditions.

The paper is organized as follows. In Section 2, we introduce our methodology and main results for statistical inference. Applications to the estimation under the spectral norm, support recovery and the estimation of latent variable graphical models are presented in Section 3. Results on linear regression are presented in Section 4 to support the main theory. Section 5 discusses possible extensions of our results and the connection between our and existing results. Numerical studies are presented in Section 6. The proof for the novel lower bound result is given in Section 7. Additional proofs are provided in Ren et al. 2015.

Methodology and statistical inference

In this section we introduce our methodology for estimating each entry and more generally, a smooth functional of any square submatrix of fixed size. Asymptotic efficiency results are stated in Section 2.3 under a sparseness assumption. The lower bound in Section 2.4 shows that the sparseness condition is sharp for the asymptotic efficiency proved in Section 2.3.

We will first introduce the methodology to estimate each entry $\omega_{ij}$, and discuss its extension to the estimation of functionals of a submatrix of the precision matrix.

The methodology is motivated by the following simple observation with $A=\{i,j\}$:

Equivalently we write a bivariate linear model

where the coefficients and error distributions are

Denote the covariance matrix of $(\eta_{i},\eta_{j})^{T}$ by

We will estimate $\Theta_{A,A}$ and expect that an efficient estimator of $\Theta_{A,A}$ yields an efficient estimation of the entries of $\Omega_{A,A}$ by inverting the estimator of $\Theta_{A,A}$.

Denote the $n$ by $p$-dimensional data matrix by $\mathbf{X.}$ The $i$th row of the data matrix is the $i$th sample $X^{(i)}$. Let $\mathbf{X}_{A}$ be the sub-matrix of $\mathbf{X}$ composed of columns indexed by $A$. Based on the regression interpretation (5), we have the following data version of the multivariate regression model

Here each row of (7) is a sample of the linear model (5). Note that $\bm{\beta}=\bm{\beta}_{A^{c},A}$ is a $p-2$ by $2$-dimensional coefficient matrix. Denote a sample version of $\Theta_{A,A}$ by

which is an oracle MLE of $\Theta_{A,A}$ based on the extra knowledge of $\bm{\beta}$. The oracle MLE of $\Omega_{A,A}$ is

Of course $\bm{\beta}$ is unknown, and we will need to estimate $\bm{\beta}$ and plug in its estimator to estimate $\bm{\varepsilon}_{A}$. This general scheme can be formally written as

where $\hat{\bm{\varepsilon}}_{A}$ is the estimated residual corresponding to a suitable estimator of $\bm{\beta}_{A^{c},A}$, that is,

Now we introduce specific estimators of $\hat{\bm{\beta}}=\hat{\bm{\beta}}_{A^{c},A}=(\hat{\bm{\beta}}_{i},\hat{\bm{\beta}}_{j})$. For each $m\in A=\{i,j\}$, we apply a scaled lasso estimator to the univariate linear regression of $\mathbf{X}_{m}$ against $\mathbf{X}_{A^{c}}$ as follows:

where $\operatorname{supp}(b)$ denotes the support of vector $b$.

Different versions of scaled lasso, in the sense of scale-free simultaneous estimation of the regression coefficients and noise level, have been considered in Städler, Bühlmann and van de Geer 2010, Antoniadis 2010 and Sun and Zhang 2010 (Sun and Zhang 2010; Sun and Zhang 2012a) among others. The $\hat{\bm{\beta}}_{m}$ in (12) is equivalent to the square-root lasso in Belloni, Chernozhukov and Wang 2011. Theoretical properties of the LSE after model selection, given in (13), were studied in Sun and Zhang 2012a (Sun and Zhang 2012a; Sun and Zhang 2013).

Our methodology can be routinely extended into a more general form. For any subset $B\subset\{1,2,\ldots,p\}$ with a bounded size, the conditional distribution of $Z_{B}$ given $Z_{B^{c}}$ is

so that the associated multivariate linear regression model is $\mathbf{X}_{B}=\penalty\mathbf{X}_{B^{c}}\bm{\beta}_{B,B^{c}}+\bm{\varepsilon}_{B}$ with $\bm{\beta}_{B^{c},B}=-\Omega_{B^{c},B}\Omega_{B,B}^{-1}$ and $\varepsilon_{B}\sim\mathcal{N}(0,\Omega_{B,B}^{-1})$. Consider a more general problem of estimating a smooth functional of $\Omega_{B,B}^{-1}$, denoted by

When $\bm{\beta}_{B^{c},B}$ is known, $\varepsilon_{B}$ is sufficient for $\Omega_{B,B}^{-1}$ due to the independence of $\varepsilon_{B}$ and $\mathbf{X}_{B^{c}}$, so that an oracle maximum likelihood estimator of $\zeta$ can be defined as

We apply an adaptive regularized estimator $\hat{\bm{\beta}}_{B^{c},B}$ by regressing $\mathbf{X}_{B}$ against $\mathbf{X}_{B^{c}}$, for example, a penalized LSE or the LSE after model selection. We estimate the residual matrix by $\hat{\bm{\varepsilon}}_{B}=\mathbf{X}_{B}-\mathbf{X}_{B^{c}}\hat{\bm{\beta}}_{B^{c},B}$, and $\zeta(\Omega_{B,B}^{-1})$ by

2 Computational complexity

For statistical inference about a single entry $\omega_{ij}$ of the precision matrix $\Omega$ with preconceived $i$ and $j$, the computational cost of the estimator (10) is of the same order as that of a single run of the scaled lasso (12).

For the estimation of the entire precision matrix $\Omega$, the definition of (10) requires the computation of ${\hat{\omega}}_{ij}$ for ${p\choose 2}$ different $A=\{i,j\}$, $i<j$. However, the computational cost for these ${p\choose 2}$ different ${\hat{\omega}}_{ij}$ is no greater than that of $(1+{\bar{s}})p$ runs of (12) where ${\bar{s}}$ is the average size of the selected model for regressing a single $\mathbf{X}_{j}$ against the other $p-1$ variables. This can be seen as follows. Define the “one-versus-rest” estimator as

and ${\hat{S}}_{i}^{(1)}=\operatorname{supp}(\hat{\bm{\beta}}_{-i,i}^{(1)})$. For $j\notin\{i\}\cup{\hat{S}}_{i}^{(1)}$, the “two-versus-rest” estimator (12) satisfies $\{\hat{\bm{\beta}}_{m},\hat{\theta}_{mm}^{1/2}\}=\{\hat{\bm{\beta}}^{(1)}_{\{i,j\}^{c},i},\sqrt{\hat{\theta}_{ii}^{(1)}}\}$ when $m=i$ and $A=\{i,j\}$. Thus we only need to carry out $1+|{\hat{S}}_{i}^{(1)}|$ runs of (12) to compute the two-versus-rest estimator $\{\hat{\bm{\beta}}_{m},\hat{\theta}_{mm}^{1/2}\}$ for all $m=i$ and $A=\{i,j\}$, $j\neq i$, where $|{\hat{S}}_{i}^{(1)}|$ denotes the cardinality of the set ${\hat{S}}_{i}^{(1)}$. Consequently, the total required runs of the scaled lasso (12) is $\sum_{i=1}^{p}(1+|{\hat{S}}_{i}^{(1)}|)=(1+{\bar{s}})p$. It follows from Theorem 11 below that $(1+{\bar{s}})p$ is of the order $\#\{(i,j)\colon\;\omega_{ij}\neq 0\}$. Thus for the computation of the estimator (10) for the entire precision matrix $\Omega$, the order of the total number of runs of (12) is the total number of edges of the graphical model corresponding to $\Omega$.

3 Statistical inference

Our analysis can be outlined as follows. We prove that estimators in the form of (10) possess the asymptotic normality and efficiency properties claimed in Theorem 1 when the following conditions hold for certain fixed constant $C_{0}$, $\varepsilon_{\Omega}\to 0$ and all $\delta\geq 1$:

with a certain complexity measure $s$ of the precision matrix $\Omega$, provided that the spectrum of $\Omega$ is bounded, and the sample size $n$ is no smaller than $(s\log p)^{2}/c_{0}$ for a sufficiently small $c_{0}>0$. This is carried out by comparing the estimator in (10) with the oracle MLE in (8) and (9) and proving

or equivalently the asymptotic normality of the oracle MLE in (9) with mean $\omega_{ij}$ and variance $n^{-1}(\omega_{ii}\omega_{jj}+\omega_{ij}^{2})$. We then prove (16), (17) and (18) for both the scaled lasso estimator (12) and the LSE after the scaled lasso selection (13). Moreover, we prove that certain thresholded versions of the proposed estimator possesses global optimality properties, as discussed below Theorem 1, under the same boundedness condition on the spectrum of $\Omega$ and a more relaxed condition on the sample size.

Suppose that conditions (16), (17) and (18) hold with $C_{0}$ and $\varepsilon_{\Omega}$. Then

with a positive constant $C_{1}$ depending on $\{C_{0},\max_{m\in A=\{i,j\}}\theta_{mm}\}$ only, and

with a constant $C_{1}^{\prime}>0$ depending on $\{c_{0}C_{1},\max_{m\in A=\{i,j\}}\{\omega_{mm},\theta_{mm}\}\}$ only.

Let $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ with $\varepsilon>0$ in (12), $\hat{\beta}_{A^{c},A}$ be the scaled lasso estimator (12) or the LSE after the scaled lasso selection (13). Then (16), (17) and (18), and thus (20) and (21), hold for all $\Omega\in\mathcal{G}^{\ast}(M,s,\lambda)$ with a certain constant $C_{0}$ depending on $\{\varepsilon,c_{0},M\}$ only and

Moreover, $Z^{\prime}$ can be defined as a linear combination of $Z_{kl}$, $kl=ii,ij,jj$.

Theorem 2 immediately yields the following results of estimation and inference for $\omega_{ij}$.

Let $\hat{\Omega}_{A,A}$ be the estimator of $\Omega_{A,A}$ in (10) with the components of $\hat{\bm{\varepsilon}}_{A}$ being the estimated residuals (11) of (12) or (13). Set $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ in (12) with certain $\delta\geq 1$ and $\varepsilon>0$. Suppose $s\leq c_{0}n/\log p$ for a sufficiently small constant $c_{0}>0$. For any small constant $\varepsilon_{0}>0$, there exists a constant $C_{2}=C_{2}(\varepsilon_{0},\varepsilon,c_{0},M)$ such that

Moreover, there exists a constant $C_{3}=C_{3}(\delta,\varepsilon,c_{0},M)$ such that

Furthermore, $\hat{\omega}_{ij}$ is asymptotically efficient with a consistent variance estimate

uniformly for all $i,j$ and $\Omega\in\mathcal{G}^{*}(M,s,\lambda)$, provided that $s=o(\sqrt{n}/\log p)$, where

The upper bounds $\max\{s\frac{\log p}{n},\sqrt{\frac{1}{n}}\}$ and $\max\{s\frac{\log p}{n},\sqrt{\frac{\log p}{n}}\}$ in equations (24) and (25), respectively, are shown to be rate-optimal in Section 2.4.

The choice of $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ is common in the literature, but can be too big and too conservative, which usually leads to some estimation bias in practice. Let $L_{n}(t)$ be the negative quantile function of $\mathcal{N}(0,1/n)$, which satisfies $L_{n}(t)\approx\sqrt{(2/n)\log p}$. In Sections 4 and 5.1 we show the value of $\lambda$ can be reduced to $(1+\varepsilon)L_{n}(k/p)$ when $\delta\vee k=o(\sqrt{n}/\log p)$.

In Theorems 2 and 3, our goal is to estimate each entry $\omega_{ij}$ of the precision matrix $\Omega$. Sometimes it is more natural to consider estimating the partial correlation $r_{ij}=-\omega_{ij}/(\omega_{ii}\omega_{jj})^{1/2}$ between $Z_{i}$ and $Z_{j}$. Let $\hat{\Omega}_{A,A}$ be estimator of $\Omega_{A,A}$ defined in (10). Our estimator of partial correlation $r_{ij}$ is defined as $\hat{r}_{ij}=-\hat{\omega}_{ij}/(\hat{\omega}_{ii}\hat{\omega}_{jj})^{1/2}$. Then the results above can be easily extended to the case of estimating $r_{ij}$. In particular, under the assumptions of Theorem 3, the estimator $\hat{r}_{ij}$ is asymptotically efficient: $\sqrt{n(1-r_{ij}^{2})^{-2}}(\hat{r}_{ij}-r_{ij})$ converges to $\mathcal{N}(0,1)$ when $s=o(\sqrt{n}/\log p)$. This asymptotic normality result was stated as Corollary $1$ in Sun and Zhang 2012b without proof.

Let $\hat{\zeta}$ be the estimator of $\zeta$ defined in (15) with the components of $\hat{\bm{\varepsilon}}_{B}$ being the estimated residuals (11) of the estimators (12) or (13). Set the penalty level $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ in (12) with certain $\delta\geq 1$ and $\varepsilon>0$. Suppose $s\leq c_{0}n/\log p$ for a sufficiently small constant $c_{0}>0$. Then

with a constant $C_{1}=C_{1}(\varepsilon,c_{0},M,|B|)$. Furthermore, $\hat{\zeta}$ is asymptotically efficient

when $\Omega\in\mathcal{G}^{*}(M,s,\lambda)$ and $s=o(\sqrt{n}/\log p)$, where $F_{\zeta}$ is the Fisher information of estimating $\zeta$ for the Gaussian model $\mathcal{N}(0,\Omega_{B,B}^{-1})$.

where $|\xi_{(1)}|\geq|\xi_{(2)}|\geq\cdots\geq|\xi_{(p)}|$. Let

Since $\xi\in B_{q}(k)$ implies $\sum_{j}\min\{1,|\xi_{j}|/\lambda\}\leq\lfloor k/\lambda^{q}\rfloor+\{q/(1-q)\}k^{1/q}\lfloor k/\penalty\lambda^{q}\rfloor^{1-1/q}/\lambda$,

when $C_{q}k_{n,p}/\lambda^{q}\leq s$, where $C_{q}=1+q2^{1/q-1}/(1-q)$ for $0<q<1$ and $C_{0}=1$. We state the extension in the following corollary.

The conclusions of Theorems 2, 3 and 4 hold with $\mathcal{G}^{*}(M,\penalty s,\lambda)$ replaced by $\mathcal{G}_{q}(M,k_{n,p})$ and $s$ by $k_{n,p}(n/\log p)^{q/2}$, $0\leq q<1$.

4 Lower bound

In this section, we derive a lower bound for estimating $\omega_{ij}$ over the matrix class $\mathcal{G}_{0}(M,k_{n,p})$ defined in (1). Assume that

for some $C_{0}>0$. Theorem 5 below implies that the assumption $k_{n,p}\frac{\log p}{n}\rightarrow 0$ is necessary for consistent estimation of any single entry of $\Omega$.

We carefully construct a finite collection of distributions $\mathcal{G}_{0}\subset\mathcal{G}_{0}(M,k_{n,p})$ and apply Le Cam’s method to show that for any estimator $\hat{\omega}_{ij}$,

for some constant $C_{1}>0$. It is relatively easy to establish the parametric lower bound $\sqrt{\frac{1}{n}}$. These two lower bounds together immediately yield Theorem 5 below.

Suppose we observe independent and identically distributed $p$-variate Gaussian random variables $X^{(1)},X^{(2)},\ldots,X^{(n)}$ with zero mean and precision matrix $\Omega=(\omega_{kl})_{p\times p}\in\mathcal{G}_{0}(M,k_{n,p})$. Under assumptions (32) and (33), we have the following minimax lower bounds:

where $c_{1},c_{2},C_{1}$, $C_{2}$, $C_{1}^{\prime}$ and $C_{2}^{\prime}$ are positive constants depending on $M$, $\nu$ and $C_{0}$ only.

The lower bound $\frac{k_{n,p}\log p}{n}$ in Theorem 5 shows that estimation of sparse precision matrix can be very different from estimation of sparse covariance matrix. The sample covariance always gives a parametric rate of estimation for every entry $\sigma_{ij}$. But for estimation of sparse precision matrix, when $k_{n,p}\gg\frac{\sqrt{n}}{\log p}$, Theorem 5 implies that it is impossible to obtain the parametric rate.

These lower bounds match the upper bounds in Corollary 1 for the proposed estimator.

Applications

where $F_{ij}=(\omega_{ii}\omega_{jj}+\omega_{ij}^{2})^{-1}$ is the Fisher information of estimating $\omega_{ij}$. The total number of edges is $p(p-1)/2$. We may apply thresholding to $\hat{\omega}_{ij}$ to correctly distinguish zero and nonzero entries. However, the variance $\omega_{ii}\omega_{jj}+\omega_{ij}^{2}$ needs to be estimated. We define the adaptive support recovery procedure as follows:

Here $\hat{\omega}_{ii}\hat{\omega}_{jj}+\hat{\omega}_{ij}^{2}$ is the natural estimate of the asymptotic variance of $\hat{\omega}_{ij}$ defined in (10), and $\xi_{0}$ is a tuning parameter which can be taken as fixed at any $\xi_{0}>2$. This thresholding estimator is adaptive. The sufficient conditions in Theorem 6 below for support recovery are much weaker than other results in literature.

Define a thresholded population precision matrix as

The following theorem shows that with high probability, ANT recovers all the strong edges without false recovery. Moreover, under the uniform signal strength condition,

Cai, Liu and Zhou 2012 showed that the rates obtained in equations (44) and (45) are optimal when $p\geq cn^{\alpha_{0}}$ for some $\alpha_{0}>1$ and $k_{n,p}=o(n^{1/2}(\log p)^{-3/2})$.

3 Estimation and inference for latent variable graphical model

Let $O$ and $H$ be two subsets of $\{1,2,\ldots,p+h\}$ with $\operatorname{Card}(O)=p$,$\operatorname{Card}(H)=h$ and $O\cup H=\{1,2,\ldots,p+h\}$. Assume that $(X_{O}^{(i)},X_{H}^{(i)})$, $i=1,\ldots,n$, are i.i.d. $(p+h)$-variate Gaussian random vectors with a positive covariance matrix $\Sigma_{(p+h)\times(p+h)}$. Denote the corresponding precision matrix by $\bar{\Omega}_{(p+h)\times(p+h)}=\Sigma_{(p+h)\times(p+h)}^{-1}$. We only have access to $\{X_{O}^{(1)},X_{O}^{(2)},\ldots,X_{O}^{(n)}\}$, while $\{X_{H}^{(1)},X_{H}^{(2)},\ldots,X_{H}^{(n)}\}$ are hidden and the number of latent components is unknown. Write $\Sigma_{(p+h)\times(p+h)}$ and $\bar{\Omega}_{(p+h)\times(p+h)}$ as follows:

where $\Sigma_{O,O}$ and $\Sigma_{H,H}$ are covariance matrices of $X_{O}^{(i)}$ and $X_{H}^{(i)}$, respectively, and from the Schur complement we have

see, for example, Horn and Johnson 1990. Define

We focus on the estimation of $\Sigma_{O,O}^{-1}$ and $S$, as the estimation of $L$ can be naturally carried out based on our results as in Chandrasekaran, Parrilo and Willsky 2012 and Ren and Zhou 2012. To make the problem identifiable we assume that $S$ is sparse, and the observed and latent variables are weakly correlated in the following sense:

which implies that both the covariance $\Sigma_{O,O}$ of observations $X_{O}^{(i)}$ and the sparse component $S=\bar{\Omega}_{O,O}$ have bounded spectrum.

With a slight abuse of notation, we denote the precision matrix $\Sigma_{O,O}^{-1}$ of $X_{O}^{(i)}$ by $\Omega$ and its inverse by $\Theta$. We propose the application of the methodology in Section 2 to i.i.d. observations $X^{(i)}$ from $\mathcal{N}(0,\Sigma_{O,O})$ with $\Omega=(s_{ij}-l_{ij})_{1\leq i,j\leq p}$ by considering the following regression:

To obtain the asymptotic normality result, condition (19) of Theorem 2 requires

with $\lambda\asymp\sqrt{(\log p)/n}$. However, when $L$ is coherent [Candès and Recht 2009] in the sense of $\{\max_{j}\sum_{i}|l_{ij}|\}^{2}\asymp p\max_{i}\sum_{j}l_{ij}^{2}\asymp p(a_{n}/n)\log p$,

Thus the conditions of Theorem 2 are not satisfied for the latent variable graphical model when $a_{n}p(\log p)^{2}\geq n$. We overcome the difficulty through a new analysis.

Let $\hat{\Omega}_{A,A}$ be the estimator of $\Omega_{A,A}$ defined in (10) with $A=\{i,j\}$ for the regression (50), where the components of $\hat{\bm{\varepsilon}}_{A}$ are the estimated residuals of (12) or (13). Let $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ for certain $\delta\geq 1$ and $\varepsilon>0$. Under assumptions (47)–(49) and $k_{n,p}\leq c_{0}n/\log p$ with a small $c_{0}$, we have

If the condition on $k_{n,p}$ is strengthened to $k_{n,p}=o(\frac{\sqrt{n}}{\log p})$, then

Let $\lambda=(1+\varepsilon)\sqrt{\frac{2\delta\log p}{n}}$ for some $\delta\geq 3$ and $\varepsilon>0$ in (12). Assume assumptions (47)–(49) hold. Then:

Under the assumptions $k_{n,p}=o(\sqrt{\frac{n}{\log p}})$ and

Regression revisited

The key element of our analysis is to establish (16), (17) and (18) for the scaled lasso estimator (12) and the LSE after the scaled lasso selection (13). The existing literature has provided theorems and arguments to carry out this task. However, several issues still require extension of existing results or explanation and modification of existing proofs. For example, the LSE after model selection is not as well understood as the lasso, and biased regression models are typically studied inexplicably, if at all. Another issue is that the penalty level used in theorems in previous sections could be too large for good numerical performance, especially for $\delta\geq 3$ in (25) of Theorems 3 and Theorems 6, 7 and 9. These issues were addressed in previous versions of this paper (\arxivurlarXiv:1309.6024) in separate lemmas. In this section, we provide a streamlined presentation of these regression results required in our analysis.

Let $\widetilde{\mathbf{X}}=(\widetilde{\mathbf{X}}_{1},\ldots,\widetilde{\mathbf{X}}_{\widetilde{p}})$ be an $n\times{\widetilde{p}}$ standardized design matrix with $\|\widetilde{\mathbf{X}}_{k}\|^{2}=n$ for all $k=1,\ldots,{\widetilde{p}}$, and $\widetilde{\mathbf{Y}}$ be a response vector satisfying

For the scaled lasso $\{\hat{\bm{\beta}}_{m},\hat{\theta}_{mm}^{1/2}\}$ in (12), $\{{\overline{\mathbf{D}}}_{A^{c}}^{1/2}\hat{\bm{\beta}}_{m},\hat{\theta}_{mm}^{1/2}\}$ can be written as

with $m\in A=\{i,j\}$, $\widetilde{\mathbf{X}}=\mathbf{X}_{A^{c}}{\overline{\mathbf{D}}}_{A^{c}}^{-1/2}$, ${\overline{\mathbf{D}}}=\operatorname{diag}(\mathbf{X}^{T}\mathbf{X}/n)$, $\widetilde{\mathbf{Y}}=\mathbf{X}_{m}$ and $\bm{\gamma}={\overline{\mathbf{D}}}_{A^{c}}^{1/2}\bm{\beta}_{m}$. For the LSE after model selection in (13), $\{{\overline{\mathbf{D}}}_{A^{c}}^{1/2}\hat{\bm{\beta}}_{m},\hat{\theta}_{mm}^{1/2}\}$ can be written as

Moreover, for both estimators, conditions (16), (17) and (18) are consequences of

To carry out an analysis of the lasso, one has to make a choice among different ways of controlling the correlations between the design and noise vectors in (53),

For $\alpha\geq 0$ and index sets $K$, the compatibility constant is defined as

We impose the following conditions on the target coefficient vector and the design:

For small penalty levels and the LSE after model selection, we also need

In (60), (61), (62) and (4), $s_{j}$ are allowed to change with $\{n,{\widetilde{p}}\}$, while $\alpha$ and $C_{j}$ are fixed constants. These conditions also make sense for deterministic designs with ${\widetilde{\varepsilon}}_{1}=0$ for deterministic conditions.

Let $k$ and $\varepsilon$ be positive real numbers and $\lambda_{0}$ be a penalty level satisfying

where $L_{n}(t)=n^{-1/2}\Phi^{-1}(1-t)$ is the $\mathcal{N}(0,1/n)$ negative quantile function. Let

We note that $L_{n}(t)=n^{-1/2}L_{1}(t)\leq\sqrt{(2/n)\log(1/t)}$ for $t\leq 1/2$, so that the right-hand side of (65) is of the order $k/\{s_{2}(\log p)^{2}\}+\sqrt{\delta/s_{2}}$. Thus condition (65) is easily satisfied even when $\varepsilon_{1}$ is a small positive number and $k$ is a moderately large number. Moreover, $\lambda$ depends on $\delta$ only through $\sqrt{\delta/s_{2}}$ in (65).

Let $\{\hat{\bm{\gamma}},\hat{\sigma}\}$ be as in (54) with data in (53) and a penalty level in (64). Let ${\widetilde{\varepsilon}}_{1}<1$ and $\lambda^{*}=L_{n-3/2}({\widetilde{\varepsilon}}_{1}/{\widetilde{p}})$. Suppose $\lambda^{*}\leq 1$ and $\delta s(\log{\widetilde{p}})/\penalty n\leq c_{0}$.

Let $s\geq s_{1}+{{s_{2}}}$, $1\leq s_{2}\leq s_{3}$, $k\geq 1$ and $\varepsilon_{1}<\varepsilon$ in (64) and (65), $\varepsilon_{1}<\varepsilon_{2}<\varepsilon$, $\alpha\geq 2(\varepsilon-\varepsilon_{2})_{+}^{-1}\{1+\varepsilon+L_{1}({\widetilde{\varepsilon}}_{1}/{\widetilde{p}})/L_{1}(k/{\widetilde{p}})\}$, $(6+e^{1/(4n-6)^{2}}){\widetilde{\varepsilon}}_{1}\leq{\widetilde{p}}^{1-\delta}{\widetilde{\varepsilon}}_{0}$ and $C_{4}\geq\sqrt{(4/L_{1}^{2}(k/{\widetilde{p}}))\log({\widetilde{p}}/{\widetilde{\varepsilon}}_{1})}/\min(\sqrt{2}-1,\varepsilon_{2}-\varepsilon_{1})$. Then there exists a constant $C_{0}$ depending on $\{\alpha,\varepsilon,\varepsilon_{2},C_{2}\}$ only such that when $C_{0}c_{0}\leq 1/2$, (60), (61), (62) and (4) imply (56), (57) and (58).

In Theorem 10, $s_{1}$ in (60) represents the complexity or the size of the coefficient vector, and $s_{2}$ represents the number of false positives we are willing to accept with the penalty level in (64). Thus $s$ is an upper bound for the total number of estimated coefficients, true or false. We summarize parallel results for the LSE after model selection as follows.

with ${\hat{S}}=\operatorname{supp}(\hat{\bm{\gamma}})$ and

Let $\lambda_{0}$ be a penalty level satisfying (64) and $\varepsilon_{1}<\varepsilon_{2}<\varepsilon_{3}<\varepsilon$. Suppose the conditions of Theorem 10 hold and that the constant factor $C_{0}$ in Theorem 10 satisfies

Then, for the parameters defined in the respective parts of Theorem 10,

If in addition, condition (61) is strengthened to

Discussion

In Theorem 2 and nearly all consequent results in Theorems 3–4 and 6–9, we have picked the penalty level $\lambda=(1+\varepsilon)\sqrt{(2\delta/n)\log p}$ for $\delta\geq 1$ ($\delta\geq 3$ for support recovery) and $\varepsilon>0$. This choice of $\lambda$ can be too conservative and may cause some finite sample estimation bias. However, in view of Theorem 10(ii) and (iii), the results in these theorems in Sections 2 and 3 still hold for penalty levels no smaller than $\lambda=(1+\varepsilon)L_{n}(k/p)\approx(1+\varepsilon)\sqrt{(2/n)\log(p/k)}$, which weakly depends on $\delta$ through (65) and the requirement of $\varepsilon>\varepsilon_{1}$.

Condition (65), with $\varepsilon<\varepsilon_{1}$, ${\widetilde{\varepsilon}}_{1}={\widetilde{p}}^{1-\delta}$ and ${\widetilde{p}}=p-2$ for the estimation of precision matrix, is the key for the choice of the smaller penalty level $\lambda=(1+\varepsilon)L_{n}(k/p)$. It provides theoretical justifications for the choice of $k\in[1,n]$ or even up to $k\asymp n\log p$ for the theory to work. Let $s_{\max}=c_{0}n/\log p$ with a sufficiently small constant $c_{0}>0$, which can be viewed as the largest possible $s\geq s_{1}+s_{2}$ in our theory. Suppose $n\leq p^{t_{0}}$ for some fixed $t_{0}<1$ and the bound $C_{3}$ for the upper sparse eigenvalue can be treated as fixed in (62) for $s_{2}\leq s_{\max}$. For $\lambda=(1+\varepsilon)\sqrt{(2/n)\log(p/k)}$ with $k\leq n\log p$ and $s_{2}\leq s_{\max}$, condition (65) can be written as

which holds for sufficiently small $ks_{\max}/(ns_{2}\log p)$. This allows $k\asymp n\log p$ for $s_{2}=s_{\max}$. For the asymptotic normality, we need $s_{2}=o(\sqrt{n}/\log p)$, so that $k=o(n^{1/2}\log p)$ is sufficient.

2 Statistical inference under unbounded condition number

3 Related works

Our methodology in this paper is related to Zhang and Zhang 2014 who proposed a LDPE approach for making inference in a high-dimensional linear model. Since $\hat{\bm{\varepsilon}}_{A}$ can be viewed as an approximate projection of $\mathbf{X}_{A}$ to the direction of $\bm{\varepsilon}_{A}$ in (10), the estimator in (10) can be viewed as an LDPE as Zhang and Zhang 2014 discussed in the regression context. See also van de Geer et al. 2014 and Javanmard and Montanari 2014. When appropriately applying their approach to our setting, their result is asymptotically equivalent to ours and also obtains the asymptotic normality. In this section, we briefly discuss their approach in the large graphical model setting.

where ${\mathbf{z}}$ is the residue after regressing $\mathbf{X}_{2}$ against all remaining columns in step one $\mathbf{X}_{A^{c}}$ using scaled lasso again. To obtain the final estimator of $\omega_{12}$, the estimator $\hat{\beta}_{2}$ of $-\omega_{12}\omega_{11}^{-1}$ should be scaled by an accurate estimator of $\omega_{11}^{-1}$, which uses the variance component of the scaled lasso estimator in the first step. It seems that two approaches are quite different. However, both approaches do the same thing: they try to estimate the partial correlation of node $Z_{1}$ and $Z_{2}$ and hence are asymptotically equivalent. Compared with their approach, our method enjoys simper form and clearer interpretation. It is worthwhile to point out that the main contribution of this paper is understanding the fundamental limit of the Gaussian graphical model in making statistical inference, which is not covered by other works.

4 Unknown mean μ\mu

In the Introduction, we assume $Z\sim\mathcal{N}(\mu,\Sigma)$ and $\mu=0$ without loss of generality. This can be seen as follows. Suppose we observe an $n\times p$ data matrix $\mathbf{X}$ with i.i.d. rows from $\mathcal{N}(\mu,\Sigma)$. Let $u^{(i)}$, $i=1,\ldots,n$, be $n$-dimensional orthonormal row vectors with $u^{(n)}=(1,\ldots,1)/\sqrt{n}$. Then $u^{(i)}\mathbf{X}$ are i.i.d. $p$-dimensional row vectors from $\mathcal{N}(0,\Sigma)$. Thus we can simply apply our methods and theory to the sample $\{u^{(i)}\mathbf{X},i=1,\ldots,n-1\}$.

Numerical studies

In this section, we present some numerical results for both asymptotic distribution and support recovery. We generate the data from $p\times p$ precision matrices with three blocks. Two cases are considered: $p=200,800$. The ratio of block sizes is $2:1:1$; that is, for a $200\times 200$ matrix, the block sizes are $100\times 100$, $50\times 50$ and $50\times 50$, respectively. The diagonal entries are $\alpha_{1},\alpha_{2},\alpha_{3}$ in three blocks, respectively, where $(\alpha_{1},\alpha_{2},\alpha_{3})=(1,2,4)$. When the entry is in the $k$th block, $\omega_{j-1,j}=\omega_{j,j-1}=0.5\alpha_{k}$, and $\omega_{j-2,j}=\omega_{j,j-2}=0.4\alpha_{k}$, $k=1,2,3$. The asymptotic variance for estimating each entry can be very different. Thus a simple procedure with a single threshold level for all entries is not likely to perform well.

We first estimate the entries in the precision matrix and partial correlations as discussed in Remark 3, and consider the distributions of these estimators. We generate a random sample of size $n=400$ from a multivariate Gaussian distribution $\mathcal{N}(0,\Sigma)$ with $\Sigma=\Omega^{-1}$. For the proposed estimators defined through (10) and (11) with the scaled lasso (12) or the LSE after model selection (13), we pick $\lambda=n^{-1/2}L_{n}(1/p)\approx\sqrt{(2/n)\log p}$; that is, $k=1$ in (64) with small adjustment in $n$ and $p$ ignored. This is justified by our theoretical results as discussed in Section 5.1.

Table 1 reports the mean and standard error of our estimators for four entries in the precision matrices and the corresponding correlations. In addition, we report the point estimates by the GLasso [Friedman, Hastie and Tibshirani 2008] and CLIME [Cai, Liu and Luo 2011] for comparison. For $p=800$, the results for the GLasso are based on 10 replications, while all other entries in the table are based on 100 replications. The GLasso is computed by the R package “glasso” with penalized diagonal (default option), while the CLIME estimators are computed by the R package “fastclime” [Pang, Liu and Vanderbei 2014]. As the GLasso and CLIME are designed for estimating precision matrices as high-dimensional objects, it is not surprising that the proposed estimator outperforms them in estimation accuracy for individual entries. Figures 1 and 2 show the histograms of the proposed estimates with the theoretical Gaussian density in Theorem 3 super-imposed. They demonstrated that the histograms match pretty well to the asymptotic distribution, especially for the LSE after model selection. The asymptotic normality leads to the following $(1-\alpha)$ confidence intervals for $\omega_{ij}$ and $r_{ij}$:

where $z_{\alpha/2}$ is the $z$-score such that $P(\mathcal{N}(0,1)>z_{\alpha/2})=\alpha/2$. Table 2 reports the empirical coverage probabilities for 95% confidence intervals, which matches well to the assigned confidence level.

Support recovery of a precision matrix is of great interest. We compare our selection results with the GLasso and CLIME. In addition to the training sample, we generate an independent sample of size 400 from the same distribution for validating the tuning parameter for the GLasso and CLIME. These estimators are computed based on the entire training sample with a range of penalty levels and a proper penalty level is chosen by minimizing the negative likelihood $\{\operatorname{trace}(\overline{\Sigma}\hat{\Omega})-\log\det(\hat{\Omega})\}$ on the validation sample, where $\overline{\Sigma}$ is the sample covariance matrix. The proposed ANT estimators are computed based on the training sample only with $\xi_{0}=2$ in the thresholding step as in (38). Tables 3 and 4 present the average selection performances as measured in the true positive, false positive and the corresponding rates. In addition to the overall performance, the summary statistics are reported for each block. The results demonstrate the selection consistency property of both ANT methods and substantial false positive for the GLasso and CLIME. It should be pointed out that the ANT takes the advantage of an additional thresholding step, while the GLasso and CLIME do not. A possible explanation of the false positive for the GLasso is a tendency for the likelihood criterion with the validation sample to pick a small penalty level. However, such an explanation seems not to hold for the CLIME, which demonstrated much lower false positive than the GLasso, as the true positive rate of the CLIME is consistently maintained at about 95% for $p=200$ and 85% for $p=800$.

Moreover, we compare the ANT with the GLasso and CLIME in a range of penalty levels. Figure 3 plots the ROC curves for the GLasso and CLIME with various penalty levels and the ANT with various thresholding levels in the follow-up procedure. It demonstrates that the CLIME outperforms the GLasso, but the two methods perform significantly more poorly than the ANT in the experiment. In addition, the circle in the plot represents the performance of the ANT with the selected threshold level as in (38). The triangle and diamond in the plot represents the performance of the GLasso and CLIME with the penalty level chosen by cross-validation, respectively. This again indicates that our method simultaneously achieves a very high true positive rate and a very low false positive rate.

Proof of Theorem 5

In this section we show that the upper bound given in Section 2.3 is indeed rate optimal. We will only establish equation (35). Equation (36) is an immediate consequence of equation (35) and the lower bound $\sqrt{\frac{\log p}{n}}$ for estimation of diagonal covariance matrices in Cai, Zhang and Zhou 2010.

Let $X^{(i)}$ be i.i.d. $\mathcal{N}(0,\Omega^{-1})$, $i=1,2,\ldots,n$, with $\Omega\in\mathcal{G}_{0}$. Let $\hat{\Omega}=(\hat{\omega}_{kl})_{p\times p}$ be an estimator of $\Omega_{m}=(\omega_{kl}^{(m)})_{p\times p}$, then

where $\alpha=\inf_{1\leq m\leq m_{\ast}}|\omega_{ij}^{(m)}-\omega_{ij}^{(0)}|$.

[Proof of Theorem 5] We shall divide the proof into three steps. Without loss of generality, consider only the cases $(i,j)=(1,1)$ and $(i,j)=(1,2)$. For the general case $\omega_{ii}$ or $\omega_{ij}$ with $i\neq j$, we could always permute the coordinates and rearrange them to the special case $\omega_{11}$ or $\omega_{12}$.

Step 1: Constructing the parameter set. We first define $\Omega_{0}$,

that is, $\Sigma_{0}=(\sigma_{kl}^{(0)})_{p\times p}$ is a matrix with all diagonal entries equal to 1, $\sigma_{12}^{(0)}=\sigma_{21}^{(0)}=b$ and the rest all zeros. Here the constant $0<b<1$ is to be determined later. For $\Omega_{m},1\leq m\leq m_{\ast}$, the construction is as follows. Without loss of generality we assume $k_{n,p}\geq 3$. Denote by $\mathcal{H}$ the collection of all $p\times p$ symmetric matrices with exactly $(k_{n,p}-2)$ elements equal to $1$ between the third and the last elements on the first row (column) and the rest all zeros. Define

where $a=\sqrt{\frac{\tau_{1}\log p}{n}}$ for some constant $\tau_{1}$ which is determined later. The cardinality of $\mathcal{G}_{0}\setminus\{\Omega_{0}\}$ is

We pick the constant $b=\frac{1}{2}(1-1/M)$ and

and prove that $\mathcal{G}_{0}\subset\mathcal{G}_{0}(M,k_{n,p})$.

For any matrix $\Omega_{m}$, $1\leq m\leq m_{\ast}$, some elementary calculations yield that

Since $b=\frac{1}{2}(1-1/M)$ and $0<\tau_{1}<\frac{(1-1/M)^{2}-b^{2}}{C_{0}}$, we have