Optimal rates of convergence for sparse covariance matrix estimation

Introduction

Minimax risk is one of the most widely used benchmarks for optimality, and substantial efforts have been made on developing minimax theories in the statistics literature. A key step in establishing a minimax theory is the derivation of minimax lower bounds and several effective lower bound arguments based on hypothesis testing have been introduced in the literature. Well-known techniques include Le Cam’s method, Assouad’s lemma and Fano’s lemma. See Le Cam (1986) and Tsybakov (2009) for more detailed discussions on minimax lower bound arguments.

Driven by a wide range of applications in high dimensional data analysis, estimation of large covariance matrices has drawn considerable recent attention. See, for example, Bickel and Levina (2008a, 2008b), El Karoui (2008), Ravikumar et al. (2008), Lam and Fan (2009), Cai, Zhang and Zhou (2010) and Cai and Liu (2011). Many theoretical results, including consistency and rates of convergence, have been obtained. However, the optimality question remains mostly open in the context of covariance matrix estimation under the spectral norm, mainly due to the technical difficulty in obtaining good minimax lower bounds.

In this paper we consider optimal estimation of sparse covariance matrices and establish the minimax rate of convergence under a range of matrix operator norm and Bregman divergence losses. A major focus is on the derivation of a rate sharp lower bound under the spectral norm loss. Conventional lower bound techniques such as the ones mentioned earlier are designed and well suited for problems with parameters that are scalar or vector-valued. They have achieved great successes in solving many nonparametric function estimation problems which can be treated exactly or approximately as estimation of a finite or infinite dimensional vector and can thus be viewed as “one-directional” in terms of the lower bound arguments. In contrast, the problem of estimating a sparse covariance matrix under the spectral norm can be regarded as a truly “two-directional” problem where one direction is along the rows and another along the columns. It cannot be essentially reduced to a problem of estimating a single or multiple vectors. As a consequence, standard lower bound techniques fail to yield good results for this matrix estimation problem. New and more general technical tools are thus needed.

In the present paper we first develop a minimax lower bound technique that is particularly well suited for treating “two-directional” problems such as estimating sparse covariance matrices. The result can be viewed as a simultaneous generalization of Le Cam’s method in one direction and Assouad’s lemma in another. This general technical tool is of independent interest and is useful for solving other matrix estimation problems such as optimal estimation of sparse precision matrices.

In the special case of $q=0$, a matrix in $\mathcal{G}_{0}(c_{n,p})$ has at most $c_{n,p}$ nonzero off-diagonal elements on each column.

The problem of estimating sparse covariance matrices under the spectral norm has been considered, for example, in El Karoui (2008), Bickel and Levina (2008b), Rothman, Levina and Zhu (2009) and Cai and Liu (2011). Thresholding methods were introduced, and rates of convergence in probability were obtained for the thresholding estimators. The parameter space $\mathcal{G}_{q}(c_{n,p})$ given in (1) also contains the uniformity class of covariance matrices considered in Bickel and Levina (2008b) as a special case. We assume that the distribution of the $X_{i}$’s is subgaussian in the sense that there is $\tau>0$ such that

Let $\mathcal{P}_{q}(\tau,c_{n,p})$ denote the set of distributions of $\mathbf{X}_{1}$ satisfying (2) and with covariance matrix $\Sigma\in\mathcal{G}_{q}(c_{n,p})$.

Our technical analysis used in establishing a rate-sharp minimax lower bound has three major steps. The first step is to reduce the original problem to a simpler estimation problem over a carefully chosen subset of the parameter space without essentially decreasing the level of difficulty. The second is to apply the general minimax lower bound technique to this simplified problem, and the final key step is to bound the total variation affinities between pairs of mixture distributions with specially designed sparse covariance matrices. The technical analysis requires ideas that are quite different from those used in the typical function/sequence estimation problems.

for $0\leq q<1$. The minimax risk of estimating the covariance matrix $\Sigma$ under the spectral norm over the class $\mathcal{P}_{q}(\tau,c_{n,p})$ satisfies

Besides the sparsity assumption considered in this paper, another commonly used structural assumption in the literature is that the covariance matrix is “bandable” where the entries decay as they move away from the diagonal. This is particularly suitable in the setting where the variables exhibit a certain ordering structure, which is often the case for time series data. Various regularization methods have been proposed and studied under this assumption. Bickel and Levina (2008a) proposed a banding estimator and obtained rate of convergence for the estimator. Cai, Zhang and Zhou (2010) established the minimax rates of convergence and introduced a rate-optimal tapering estimator. In particular, Cai, Zhang and Zhou (2010) derived rate sharp minimax lower bounds for estimating bandable matrices. It should be noted that the lower bound techniques used there do not lead to a good result for estimating sparse covariance matrices under the spectral norm.

General lower bound for minimax risk

In this section we develop a new general minimax lower bound technique that is particularly well suited for treating “two-directional” problems such as estimating sparse covariance matrices. The new method can be viewed as a generalization of both Le Cam’s method and Assouad’s lemma. To help motivate and understand the new lower bound argument, it is useful to briefly review Le Cam’s method and Assouad’s lemma.

Write $\Theta_{1}=\{\theta_{1},\ldots,\theta_{D}\}$. One can view the lower bound in (5) as obtained from testing the simple hypothesis $H_{0}\dvtx\theta=\theta_{0}$ against the composite alternative $H_{1}\dvtx\theta\in\Theta_{1}$.

Assouad’s lemma works with a hypercube $\Theta=\{0,1\}^{r}$. It is based on testing a number of pairs of simple hypotheses and is connected to multiple comparisons. For a parameter $\theta=(\theta_{1},\ldots,\theta_{r})$ where $\theta_{i}\in\{0,1\}$, one tests whether $\theta_{i}=0$ or $1$ for each $1\leq i\leq r$ based on the observation $X$. For each pair of simple hypotheses, there is a certain loss for making an error in the comparison. The lower bound given by Assouad’s lemma is a combination of losses from testing all pairs of simple hypotheses. Let

be the Hamming distance on $\Theta$. Assouad’s lemma gives a lower bound for the maximum risk over the hypercube $\Theta$ of estimating an arbitrary quantity $\psi(\theta)$ belonging to a metric space with metric $d$.

In comparison, the standard lower bound arguments work with either $\Gamma$ or $\Lambda$ alone. For example, Assouad’s lemma considers only the parameter set $\Gamma$ and the Le Cam’s method typically applies to a parameter set like $\Lambda$ with $r=1$. For $\theta=(\gamma,\lambda)\in\Theta$, denote the projection of $\theta$ to $\Gamma$ by $\gamma(\theta)=\gamma$ and to $\Lambda$ by $\lambda(\theta)=\lambda$.

The following lemma gives a lower bound for the maximum risk over the parameter set $\Theta$ of estimating a functional $\psi(\theta)$ belonging to a metric space with metric $d$.

In applications of Lemma 3, for a $\gamma=(\gamma_{1},\ldots,\gamma_{r})\in\Gamma$ where $\gamma_{i}$ takes value or $1$, and a $\lambda=(\lambda_{1},\ldots,\lambda_{r})\in\Lambda$ where each $\lambda_{i}\in B$ is a $p$-dimensional nonzero row vector, the element $\theta=(\gamma,\lambda)\in\Theta$ can be equivalently viewed as an $r\times p$ matrix

Note that the lower bound (10) reduces to the classical Assouad lemma when $\Lambda$ contains only one matrix for which every row is nonzero, and becomes a two-point argument of Le Cam with one point against a mixture when $r=1$. The proof of this lemma is given in Section 7. The technical argument is an extension of that of Assouad’s lemma. See Assouad (1983), Yu (1997) and van der Vaart (1998).

The advantage of this method is the ability to break down the lower bound calculations for the whole matrix estimation problem into calculations for individual rows so that the overall analysis is simplified and more tractable. Although the tool is introduced here for the purpose of estimating a sparse covariance matrix, it is of independent interest and is expected to be useful for solving other matrix estimation problems as well.

Bounding the total variation affinity between two mixture distributions in (10) is quite challenging in general. The following well-known result on the affinity is helpful in some applications. It provides lower bounds for the affinity between two mixture distributions in terms of the affinities between simpler distributions in the mixtures.

More specifically, in our construction of the parameter set for establishing the minimax lower bound, $r$ is the number of possibly nonzero rows in the upper triangle of the covariance matrix, and $\Lambda$ is the set of matrices with $r$ rows to determine the upper triangle matrix. Recall that the projection of $\theta\in\Theta$ to $\Gamma$ is $\gamma(\theta)=\gamma=(\gamma_{i}(\theta))_{1\leq i\leq r}$ and the projection of $\theta$ to $\Lambda$ is $\lambda(\theta)=\lambda=(\lambda_{i}(\theta))_{1\leq i\leq r}$. More generally, for a subset $A\subseteq\{1,2,\ldots,r\}$, we define a projection of $\theta$ to a subset of $\Gamma$ by $\gamma_{A}(\theta)=(\gamma_{i}(\theta))_{i\in A}$. A particularly useful example of set $A$ is

for which $\gamma_{\{-i\}}(\theta)=(\gamma_{1}(\theta),\ldots,\gamma_{i-1}(\theta),\gamma_{i+1}(\theta),\gamma_{r}(\theta))$ and in this case for convenience we set $\gamma_{-i}=\gamma_{\{-i\}}$. $\lambda_{A}(\theta)$ and $\lambda_{-i}(\theta)$ are defined similarly. We also define the set $\Lambda_{A}=\{\lambda_{A}(\theta)\dvtx\theta\in\Theta\}$. A special case is $A=\{-i\}$.

Now we define a subset of $\Theta$ to reduce the problem of estimating $\Theta$ to a problem of estimating $\lambda_{i}(\theta)$. For $a\in\{0,1\}$, $b\in\{0,1\}^{r-1}$ and $c\in\Lambda_{-i}\subseteq B^{r-1}$, let

and $D_{(i,b,c)}=\operatorname{Card}(\Theta_{(i,a,b,c)})$. Note that the cardinality of $\Theta_{(i,a,b,c)}$ on the right-hand side does not depend on the value of $a$ due to the Cartesian product structure of $\Theta=\Gamma\otimes\Lambda$. Define the mixture distribution

The parameter $\theta$ is seen uniformly distributed over $\Theta$. Let

and an average of $h(\gamma_{-i},\lambda_{-i})$ over the set $\Theta_{-i}$ is defined as follows:

where the distribution of $(\gamma_{-i},\lambda_{-i})$ is induced by the uniform distribution over $\Theta$.

Lower bound for estimating sparse covariance matrix under the spectral norm

We now state and prove the minimax lower bound for estimating a sparse covariance matrix over the parameter space $\mathcal{G}_{q}(c_{n,p})$ under the spectral norm. The derivation of the lower bounds relies heavily on the general lower bound technique developed in the previous section. It also requires a careful construction of a finite subset of the parameter space and detailed calculations of an effective lower bound for the total variation affinities between mixtures of multivariate Gaussian distributions.

Let $\mathbf{X}_{1},\ldots,\mathbf{X}_{n}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(\mu,\Sigma_{p\times p})$. The minimax risk for estimating the covariance matrix $\Sigma$ over the parameter space $\mathcal{G}_{q}(c_{n,p})$ with $c_{n,p}\leq Mn^{{(1-q)}/{2}}(\log p)^{-{(3-q)}/{2}}$ satisfies

for some constant $c>0$, where $|\!|\!|\cdot|\!|\!|$ denotes the matrix spectral norm.

Theorem 2 yields immediately a minimax lower bound for the more general subgaussian case under assumption (2),

It has been shown in Cai, Zhang and Zhou (2010) that

by constructing a parameter space with only diagonal matrices. It then suffices to show that

The proof of Theorem 2 contains three major steps. In the first step we construct in detail a finite subset $\mathcal{F}_{*}$ of the parameter space $\mathcal{G}_{q}(c_{n,p})$ such that the difficulty of estimation over $\mathcal{F}_{*}$ is essentially the same as that of estimation over $\mathcal{G}_{q}(c_{n,p})$. The second step is the application of Lemma 3 to the carefully constructed parameter set $\mathcal{F}_{*}$. Finally in the third step we calculate the factor $\alpha$ defined in (11) and the total variation affinity between two multivariate normal mixtures. Bounding the affinity is technically involved. The main ideas of the proof are outlined here, and detailed proofs of some technical lemmas used here are deferred to Section 7.

Proof of Theorem 2 The proof is divided into three main steps.

Step 1: Constructing the parameter set. Let $r=\lfloor p/2\rfloor$, where $\lfloor x\rfloor$ denotes the largest integer less than or equal to $x$, and let $B$ be the collection of all row vectors $b=(v_{j})_{1\leq j\leq p}$ such that $v_{j}=0$ for $1\leq j\leq p-r$ and $v_{j}=0$ or $1$ for $p-r+1\leq j\leq p$ under the constraint the total number of 1s is $\|b\|_{0}=k$, where the value of $k$ will be specified later. We shall treat each $(b_{1},\ldots,b_{r})\in B^{r}$ as an $r\times p$ matrix with the $i$th row equal to $b_{i}$.

Set $\Gamma=\{0,1\}^{r}$. Define $\Lambda\subset B^{r}$ to be the set of all elements in $B^{r}$ such that each column sum is less than or equal to $2k$. For each component $\lambda_{m}$, $1\leq m\leq r$, of $\lambda=(\lambda_{1},\ldots,\lambda_{r})\in\Lambda$, define a $p\times p$ symmetric matrix $A_{m}(\lambda_{m})$ by making the $m$th row of $A_{m}(\lambda_{m})$ equal to $\lambda_{m}$, the $m$th column equal to $\lambda_{m}^{T}$ and the rest of the entries . Note that for each $\lambda=(\lambda_{1},\ldots,\lambda_{r})\in\Lambda$, each column/row sum of the matrix $\sum_{m=1}^{r}A_{m}(\lambda_{m})$ is less than or equal to $2k$.

It is easy to see that in the Gaussian case $|\!|\!|\Sigma_{p\times p}|\!|\!|\leq\tau$ is a sufficient condition for (2). Without loss of generality we assume that $\tau>1$ in the subgaussianity assumption (2); otherwise we replace $I_{p}$ in (19) by $cI_{p}$ with a small constant $c>0$. Finally we define a collection $\mathcal{F}_{\ast}$ of covariance matrices as

Note that each $\Sigma\in\mathcal{F}_{\ast}$ has value $1$ along the main diagonal, and contains an $r\times r$ submatrix, say, $A$, at the upper right corner, $A^{T}$ at the lower left corner and elsewhere. Each row of $A$ is either identically (if the corresponding $\gamma$ value is ) or has exactly $k$ nonzero elements with value $\epsilon_{n,p}$.

We now specify the values of $\epsilon_{n,p}$ and $k$ to ensure $\mathcal{F}_{\ast}\subset\mathcal{G}_{q}(c_{n,p})$. Set $\epsilon_{n,p}=\upsilon\sqrt{\frac{\log p}{n}}$ for a fixed small constant $\upsilon$, and let $k=\max(\lceil\frac{1}{2}c_{n,p}\epsilon_{n,p}^{-q}\rceil-1,0)$ which implies

Note that $\epsilon_{n,p}$ and $k$ satisfy

and consequently every $\Sigma(\theta)$ is diagonally dominant and positive definite, and $|\!|\!|\Sigma(\theta)|\!|\!|\leq|\!|\!|\Sigma(\theta)|\!|\!|_{1}\leq 2k\epsilon_{n,p}+1<\tau$. Thus we have $\mathcal{F}_{\ast}\subset\mathcal{G}_{q}(c_{n,p})$, and the subgaussianity assumption (2) is satisfied.

For $\alpha$ defined in equation (24) we have

The key technical difficulty is in bounding the affinity between the Gaussian mixtures. The proof is quite involved.

Finally, the minimax lower bound for estimation over $\mathcal{G}_{q}(c_{n,p})$ is obtained by putting together the bounds given in Lemmas 5 and 6,

Minimax upper bound under the spectral norm

Section 3 developed a minimax lower bound for estimating a sparse covariance matrix under the spectral norm over $\mathcal{G}_{q}(c_{n,p})$. In this section we shall show that the lower bound is rate-sharp and therefore establish the optimal rate of convergence. To derive a minimax upper bound, we shall consider the properties of a thresholding estimator introduced in Bickel and Levina (2008b). Given a random sample $\{\mathbf{X}_{1},\ldots,\mathbf{X}_{n}\}$ of $p$-variate observations drawn from a distribution in $\mathcal{P}_{q}(\tau,c_{n,p})$, the sample covariance matrix is

which is an unbiased estimate of $\Sigma$, and the maximum likelihood estimator of $\Sigma$ is

when $\mathbf{X}_{l}$’s are normally distributed. These two estimators are close to each other for large $n$. We shall construct estimators of the covariance matrix $\Sigma$ by thresholding the maximum likelihood estimator $\Sigma^{\ast}$.

Note that the subgaussianity condition (2) implies

Then the empirical covariance $\sigma_{i,j}^{\ast}$ satisfies the following large deviation result that there exist constants $C_{1}>0$ and $\gamma>0$ such that

for $|t|\leq\delta$, where $C_{1},$ $\gamma$ and $\delta$ are constants and depend only on $\tau$. See Saulis and Statulevičius (1991) and Bickel and Levina (2008a). Inequality (26) implies $\sigma_{ij}^{\ast}$ behaves like a subgaussian random variable. In particular for $t=\gamma\sqrt{\frac{\log p}{n}}$ we have

Define the thresholding estimator $\hat{\Sigma}=(\hat{\sigma}_{ij})_{p\times p}$ by

This thresholding estimator was first proposed in Bickel and Levina (2008b) in which a rate of convergence of the loss function in probability was given over the uniformity class $\mathcal{G}_{q}^{\ast}(c_{n,p})$. Here we provide an upper bound for mean squared spectral norm error over the parameter space $\mathcal{G}_{q}(c_{n,p})$.

Throughout the rest of the paper we denote by $C$ a generic positive constant which may vary from place to place. The following theorem shows that the thresholding estimator defined in (28) is rate optimal over the parameter space $\mathcal{G}_{q}(c_{n,p})$.

The thresholding estimator $\hat{\Sigma}$ given in (28) satisfies, for some constant $C>0$,

Consequently, the minimax risk of estimating the sparse covariance matrix $\Sigma$ over $\mathcal{G}_{q}(c_{n,p})$ satisfies

Theorem 3 shows that the optimal rate of convergence for estimating a sparse covariance matrix over $\mathcal{G}_{q}(c_{n,p})$ under the squared spectral norm is $c_{n,p}^{2}(\frac{\log p}{n})^{1-q}$. In Bickel and Levina (2008b) the uniformity class $\mathcal{G}_{q}^{\ast}(c_{n,p})$ defined in (16) was considered. We shall now show that the same minimax rate of convergence holds for estimation over $\mathcal{G}_{q}^{\ast}(c_{n,p})$. It is easy to check in the proof of the lower bound that for every $\Sigma\in\mathcal{F}_{\ast}$ defined in (20), we have

The minimax risk for estimating the covariance matrix under the spectral norm over the uniformity class $\mathcal{G}_{q}^{\ast}(c_{n,p})$ satisfies

The thresholding estimator $\hat{\Sigma}$ defined by (28) is positive definite with high probability, but it is not guaranteed to be positive definite. A simple additional step can make the final estimator positive semi-definite and achieve the optimal rate of convergence. Write the eigen-decomposition of $\hat{\Sigma}$ as

where $\hat{\lambda}_{i}$’s and $v_{i}$’s are the eigenvalues and eigenvectors of $\hat{\Sigma}$, respectively. Let $\hat{\lambda}_{i}^{+}=\max(\hat{\lambda}_{i},0)$ be the positive part of $\hat{\lambda}_{i}$ and define

The resulting estimator $\hat{\Sigma}^{+}$ is positive semi-definite and attains the same rate as the original thresholding estimator $\hat{\Sigma}$. This method can be applied to the tapering estimator in Cai, Zhang and Zhou (2010) as well to make the estimator positive semi-definite, while still achieving the optimal rate.

Optimal estimation under Bregman divergences

We have so far focused on the optimal rate of convergence under the spectral norm. In this section we turn to minimax estimation of sparse covariance matrices under a class of Bregman divergence losses which include Stein’s loss, Frobenius norm and von Neumann’s entropy as special cases. Bregman matrix divergences have been used for matrix estimation and matrix approximation problems; see, for example, Dhillon and Tropp (2007), Ravikumar et al. (2008) and Kulis, Sustik and Dhillon (2009). In this section we establish the optimal rate of convergence uniformly for a class of Bregman divergence losses.

Bregman (1967) introduced the Bregman divergence as a dissimilarity measure between vectors,

where $X$ and $Y$ are real symmetric matrices, and $\phi$ is a differentiable strictly convex function over the space. See Censor and Zenios (1997) and Kulis, Sustik and Dhillon (2009). A particularly interesting class of $\phi$ is

$\varphi(\lambda)=$ $-\log\lambda,$ or equivalently $\phi(X)=-\log\det(X)$. The corresponding Bregman divergence can be written as

which is often called Stein’s loss in the statistical literature.

$\varphi(\lambda)=$ $\lambda\log\lambda-\lambda,$ or equivalently $\phi(X)=\operatorname{tr}(X\log X-X)$, where $X$ is positive definite such that $\log X$ is well defined. The corresponding Bregman divergence is the von Neumann divergence

$\varphi(\lambda)=$ $\lambda^{2},$ or equivalently $\phi(X)=\operatorname{tr}(X^{2})$. The resulting Bregman divergence is the squared Frobenius norm

for $X=(x_{ij})_{1\leq i,j\leq p}$ and $Y=(y_{ij})_{1\leq i,j\leq p}$.

Define a class $\Psi$ of functions $\varphi$ satisfying the following conditions:

$\varphi$ is twice differentiable, real-valued and strictly convex over $\lambda\in(0,\infty)$;

$|\varphi(\lambda)|\leq C\lambda^{r}$ for some $C>0$ and some real number $r$ uniformly over $\lambda\in(0,\infty)$;

For every positive constants $\epsilon_{2}$ and $M_{2}$ there are some positive constants $c_{L}$ and $c_{u}$ depending on $\epsilon_{2}$ and $M_{2}$ such that $c_{L}\leq\varphi^{{}^{\prime\prime}}(\lambda)\leq c_{u}$ for all $\lambda\in[\epsilon_{2},M_{2}]$.

In this paper, we shall consider the following class of Bregman divergences:

It is easy to see that Stein’s loss, von Neumann’s divergence and the squared Frobenius norm are in this class.

Let $\epsilon_{1}>0$ be a positive constant. Let $\mathcal{P}_{q}^{B}(\tau,c_{n,p})$ denote the set of distributions of $\mathbf{X}_{1}$ satisfying (2) and with covariance matrix

The following theorem gives a unified result on the minimax rate of convergence for estimating the covariance matrix over the parameter space $\mathcal{P}_{q}^{B}(\tau,c_{n,p})$ for all Bregman divergences $\phi\in\Phi$ defined in (32).

Assume that $c_{n,p}\leq Mn^{{(1-q)}/{2}}(\log p)^{-{{(3-q)}/{2}}}$ for some $M>0$ and $0\leq q<1$. The minimax risk over $\mathcal{P}_{q}^{B}(\tau,c_{n,p})$ under the loss function

for all Bregman divergences $\phi\in\Phi$ defined in (32) satisfies

Note that Theorem 4 gives the minimax rate of convergence uniformly under all Bregman divergences defined in (32). For an individual Bregman divergence loss, the condition that all eigenvalues are bounded away from is not needed if the function $\varphi$ is well behaved at . For example, such is the case for the Frobenius norm.

The optimal rate of convergence is attained by a modified thresholding estimator. Let $\hat{\Sigma}=(\hat{\sigma}_{ij})_{1\leq i,j\leq p}$ be the thresholding estimator given in (28). Define the final estimator of $\Sigma$ by

Let $\mathcal{P}_{q}^{\ast B}(\tau,c_{n,p})$ denote the set of distributions of $\mathbf{X}_{1}$ satisfying (2) and with covariance matrix $\Sigma\in\mathcal{G}_{q}^{\ast B}(c_{n,p})=\mathcal{G}_{q}^{\ast}(c_{n,p})\cap\{\Sigma\dvtx\lambda_{min}\geq\epsilon_{1}\}$. Then under the same conditions as in Theorem 4,

Discussions

Moreover, the thresholding estimator $\hat{\Sigma}$ defined in (28) is rate-optimal.

The spectral norm of a matrix depends on the entries in a subtle way and the “interactions” among different rows/columns must be taken into account. The lower bound argument developed in this paper is aimed at treating “two-directional” problems by mixing over both rows and columns. It can be viewed as a simultaneous application of Le Cam’s method in one direction and Assouad’s lemma in another. In contrast, for sequence estimation problems, we typically need one or the other, but not both at the same time. The lower bound techniques developed in this paper can be used to solve other matrix estimation problems. For example, Cai, Liu and Zhou (2011) applied the general lower bound argument to the problem of estimating sparse precision matrices under the spectral norm and established the optimal rate of convergence. This problem is closely connected to graphical model selection. The derivations of both the lower and upper bounds are involved. For reasons of space, we shall report the results elsewhere.

In addition to the hard thresholding estimator used in Bickel and Levina (2008b), Rothman, Levina and Zhu (2009) considered a class of thresholding rules with more general thresholding functions, including soft thresholding and adaptive Lasso. It is straightforward to show that these thresholding estimators with the same choice of threshold level used in (28) also attains the optimal rate of convergence over the parameter space $\mathcal{G}_{q}(c_{n,p})$ under mean squared spectral norm error as well as under the class of Bregman divergence losses considered in Section 5 with the same modification as in (34). Therefore, the choice of the thresholding function is not important as far as the rate optimality is concerned.

Proofs

In this section we prove the general lower bound result given in Lemma 3, Theorems 3 and 4 as well as some of the important technical lemmas used in the proof of Theorem 2 given in Section 3. The proofs of a few technical results used in this section are deferred to the supplementary material [Cai and Zhou (2012)]. Throughout this section, we denote by $C$ a generic constant that may vary from place to place.

We first bound the maximum risk by the average over the whole parameter set,

Set $\hat{\theta}=\arg\min_{\theta\in\Theta}d^{s}(T,\psi(\theta))$. Note that the minimum is not necessarily unique. When it is not unique, pick $\hat{\theta}$ to be any point in the minimum set. Then the triangle inequality for the metric $d$ gives

where the last inequality is due to the fact $d(\psi(\hat{\theta}),T)=d(T,\psi(\hat{\theta}))\leq d(T,\psi(\theta))$ from the definition of $\hat{\theta}$. Equations (7.1) and (7.1) together yield

where the last step follows from the definition of $\alpha$ in equation (11).

The right-hand side can be further written as

The following elementary result is useful to establish the lower bound for the minimax risk. See, for example, page 40 of Le Cam (1973).

2 Proof of Lemma 5

Let $v=(v_{i})$ be a column $p$-vector with $v_{i}=0$ for $1\leq i\leq p-r$ and $v_{i}=1$ for $p-r+1\leq i\leq p$, that is, $v=(1\{p-r+1\leq i\leq p\})_{p\times 1}$. Set $w=(w_{i})=[\Sigma(\theta)-\Sigma(\theta^{\prime})]v$. Note that for each $i$, if $|\gamma_{i}(\theta)-\gamma_{i}(\theta^{\prime})|=1$, we have $|w_{i}|=k\epsilon_{n,p}$. Then there are at least $H(\gamma(\theta),\gamma(\theta^{\prime}))$ number of elements $w_{i}$ with $|w_{i}|=k\epsilon_{n,p}$, which implies

Since $\|v\|^{2}=r\leq p$, the equation above yields

when $H(\gamma(\theta),\gamma(\theta^{\prime}))\geq 1$.

3 Proof of Lemma 6

(i) There exists a constant $c_{2}<1$ such that

Here $\mathbf{S}_{(p-1)\times(p-1)}=(s_{ij})_{2\leq i,j\leq p}$ is a symmetric matrix uniquely determined by $(\gamma_{-1},\lambda_{-1})=((\gamma_{2},\ldots,\gamma_{r}),(\lambda_{2},\ldots,\lambda_{r}))$ where for $i\leq j$,

where $\|\mathbf{r}\|_{0}=k$ with nonzero elements of $r$ equal $\epsilon_{n,p}$ and the submatrix $\mathbf{S}_{(p-1)\times(p-1)}$ is the same as the one for $\Sigma_{0}$ given in (41).

The following lemma is useful for calculating the cross product terms in the chi-squared distance between Gaussian mixtures. The proof of the lemma is straightforward and is thus omitted.

Let $g_{i}$ be the density function of $N(0,\Sigma_{i})$ for $i=0,1$ and $2$, respectively. Then

Let $\Sigma_{0}$ be defined in (41) and determined by $(\gamma_{-1},\lambda_{-1})$. Let $\Sigma_{1}$ and $\Sigma_{2}$ be of the form (42) with the first row $\lambda_{1}$ and $\lambda_{1}^{\prime}$, respectively. Set

We sometimes drop the indices $(\lambda_{1}$, $\lambda_{1}^{\prime})$ and $(\gamma_{-1},\lambda_{-1})$ from $\Sigma_{i}$ to simplify the notation whenever there is no ambiguity. Then each term in the chi-squared distance on the left-hand side of (40) can be expressed as in the form of

It is a subset of $\Theta_{-1}$ in which the element can pick both $a_{1}$ and $a_{2}$ as the first row to form parameters in $\Theta$. From Lemma 9 the average of the chi-squared distance on the left-hand side of equation (40) can now be written as

where $\lambda_{1}$ and $\lambda_{1}^{\prime}$ are independent and uniformly distributed over $\Lambda_{1}(\lambda_{-1})$ (not over $B$) for given $\lambda_{-1}$, and the distribution of $(\gamma_{-1},\lambda_{-1})$ given $(\lambda_{1},\lambda_{1}^{\prime})$ is uniform over $\Theta_{-1}$ $(\lambda_{1},\lambda_{1}^{\prime})$, but the marginal distribution of $\lambda_{1}$ and $\lambda_{1}^{\prime}$ are not independent and uniformly distributed over $B$.

Let $\Sigma_{1}$ and $\Sigma_{2}$ be two covariance matrices of the form (42). Note that $\Sigma_{1}$ and $\Sigma_{2}$ differ from each other only in the first row/column. Then $\Sigma_{i}-\Sigma_{0}$, $i=1$ or $2$, has a very simple structure. The nonzero elements only appear in the first row/column, and in total there are at most $2k$ nonzero elements. This property immediately implies the following lemma which makes the problem of studying the determinant in Lemma 9 relatively easy. The proof of Lemma 10 below is given in the supplementary material.

Let $\Sigma_{0}$ be defined in (41) and let $\Sigma_{1}$ and $\Sigma_{2}$ be two covariance matrices of the form (42). Define $J$ to be the number of overlapping $\epsilon_{n,p}$’s between $\Sigma_{1}$ and $\Sigma_{2}$ on the first row, and

There are index subsets $I_{r}$ and $I_{c}$ in $\{2,\ldots,p\}$ with $\operatorname{Card}(I_{r})=\operatorname{Card}(I_{c})=k$ and $\operatorname{Card}(I_{r}\cap I_{c})=J$ such that

and the matrix $(\Sigma_{0}-\Sigma_{1})(\Sigma_{0}-\Sigma_{2})$ has rank $2$ with two identical nonzero eigenvalues $J\epsilon_{n,p}^{2}$ when $J>0$.

The matrix $Q$ is determined by two interesting parts, the first element $q_{11}=J\epsilon_{n,p}^{2}$ and a very special $k\times k$ square matrix $(q_{ij}\dvtx i\in I_{r}\mbox{ and }j\in I_{c})$ with all elements equal to $\epsilon_{n,p}^{2}$. The following result, which is proved in the supplementary material, shows that $R_{\lambda_{1},\lambda_{1}^{\prime}}^{\gamma_{-1},\lambda_{-1}}$ is approximately equal to

Let $R_{\lambda_{1},\lambda_{1}^{\prime}}^{\gamma_{-1},\lambda_{-1}}$ be defined in equation (43). Then

where $R_{1,\lambda_{1},\lambda_{1}^{\prime}}^{\gamma_{-1},\lambda_{-1}}$ satisfies, uniformly over all $J$,

With the preparations given above, we are now ready to establish equation (40) and thus complete the proof of Lemma 6.

Recall that $J$ is the number of overlapping $\epsilon_{n,p}$’s between $\Sigma_{1}$ and $\Sigma_{2}$ on the first row. It is easy to see that $J$ has the hypergeometric distribution as $\lambda_{1}$ and $\lambda_{1}^{\prime}$ vary in $B$ for each given $\lambda_{-1}$. For $0\leq j\leq k$,

where $\frac{k!}{(k-j)!}$ is a product of $j$ term with each term $\leq k$ and for $\frac{p_{\lambda_{-1}}!(p_{\lambda_{-1}}-2k+j)!}{[(p_{\lambda_{-1}}-k)!]^{2}}$ it is bounded below by a product of $j$ term with each term $\geq p_{\lambda_{-1}}-j$. Since $p_{\lambda_{-1}}\geq p/4-1$ for all $\lambda_{-1}$, we have

by setting $c_{2}^{2}=3/4$, where the last step follows from $\upsilon^{2}<\frac{\beta-1}{54\beta}$ and $k^{2}=O(\frac{n}{\log p})=O(\frac{p^{1/\beta}}{\log p})$ as defined in Section 3.

The condition $p\geq n^{\beta}$ for some $\beta>1$ is assumed so that

for some $\varepsilon>0$ to make the term (7.3) to be $o(1)$.

4 Proof of Theorem 3

The following lemma, which is proved in Cai and Zhou (2009), is now useful to prove Theorem 3.

Let $D=(d_{ij})_{1\leq i,j\leq p}$ with $d_{ij}=(\hat{\sigma}_{ij}-\sigma_{ij})I(A_{ij}^{c})$. Then

Set $k^{\ast}=\lfloor c_{n,p}({\frac{n}{\log p}})^{q/2}\rfloor$. Then we have

Putting $R_{1}$ and $R_{2}$ together yields that for some constant $C>0$ ,

Theorem 3 is proved by combining equations (7.4), (7.4) and (52).

5 Proof of Theorem 4

We establish separately the lower and upper bounds under the Bregman divergence losses. The following lemma relates a general Bregman divergence to the squared Frobenius norm.

Assume that all eigenvalues of two symmetric matrices $X$ and $Y$ belong to $[\epsilon_{2},M_{2}]$. Then there exist constants $c_{2}>c_{1}>0$ depending only on $\epsilon_{2}$ and $M_{2}$ such that for all $\phi\in\Phi$ defined in (32),

Let the eigen decompositions of $X$ and $Y$ be

For every $\phi(X)=\sum_{i=1}^{p}\varphi(\lambda_{i})$ it is easy to see that

See Kulis, Sustik and Dhillon (2009), Lemma 1. The Taylor expansion gives

where $\xi_{ij}$ is in between $\lambda_{i}$ and $\gamma_{j}$ and then contained in $[\epsilon_{2},M_{2}]$. From the assumption in (32), there are constants $c_{L}$ and $c_{u}$ such that $c_{L}\leq\varphi^{\prime\prime}(\lambda)\leq c_{u}$ for all $\lambda$ in $[\epsilon_{2},M_{2}]$, which immediately implies

Lower bound under Bregman matrix divergences. It is trivial to see that

by constructing a parameter space with only diagonal matrices. It is then enough to show that there exists some constant $c>0$ such that

for all $\phi\in\Phi$ defined in (32). Equation (53) implies

Convexity of $\varphi$ implies $\varphi(\lambda_{i})-\varphi(\gamma_{j})-\varphi^{\prime}(\gamma_{j})\cdot(\lambda_{i}-\gamma_{j})$ is nonnegative and increasing when $\lambda_{i}$ moves away from the range $[\epsilon_{1},2\tau]$ of those eigenvalues $\gamma_{j}$’s of $\Sigma(\theta)$. From Lemma 13 there is a universal constant $c_{L}$ such that

where the last equality is from the same argument for equation (54).

It then suffices to study the lower bound under the Frobenius norm. Similar to the lower bound under the spectral norm one has

and it follows from Lemma 6 that there is a constant $c>0$ such that

Upper bound under Bregman matrix divergences. We now show that there exists an estimator $\hat{\Sigma}$ such that

some constant $c>0$, uniformly over all $\phi\in\Phi$ and $\Sigma\in\mathcal{P}_{q}^{B}(\tau,c_{n,p})$. Let $A_{0}=\bigcap_{i,j}A_{ij}$, where $A_{ij}$ is defined in (49). Lemma 12 yields that

Let $\hat{\Sigma}_{B}$ be defined in equation (34). Then for all $\Sigma\in\mathcal{G}_{q}^{B}(\rho,c_{n,p})$

Write $\hat{\Sigma}_{B}=\Sigma+(\hat{\Sigma}_{B}-\Sigma)$. Since $|\!|\!|\hat{\Sigma}_{B}-\Sigma|\!|\!|\leq|\!|\!|\hat{\Sigma}_{B}-\Sigma|\!|\!|_{1}$, the lemma is then a direct consequence of Lemma 12 and equation (7.4) which implies $|\!|\!|\hat{\Sigma}_{B}-\Sigma|\!|\!|_{1}\leq Cc_{n,p}(\frac{\log p}{n})^{(1-q)/2}\rightarrow 0$ over $A_{0}$.

The second term in (7.5) is negligible since

by applying the Cauchy–Schwarz inequality twice. We now consider the first term in equation (7.5). Set $k^{\ast}=\lfloor c_{n,p}({\frac{n}{\log p}})^{q/2}\rfloor$. Then we have

Supplement to “Optimal rates of convergence for sparse covariance matrix estimation” \slink[doi]10.1214/12-AOS998SUPP \sdatatype.pdf \sfilenameaos998_supp.pdf \sdescriptionIn this supplement we prove the additional technical lemmas used in the proof of Lemma 6.

References