Minimax sparse principal subspace estimation in high dimensions

Introduction

Principal components analysis (PCA) was introduced in the early 20th century [Pearson (1901), Hotelling (1933)] and is arguably the most well known and widely used technique for dimension reduction. It is part of the mainstream statistical repertoire and is routinely used in numerous and diverse areas of application. However, contemporary applications often involve much higher-dimensional data than envisioned by the early developers of PCA. In such high-dimensional situations, where the number of variables $p$ is of the same order or much larger than the number of observations $n$, serious difficulties emerge: standard PCA can produce inconsistent estimates of the principal directions of variation and lead to unreliable conclusions [Johnstone and Lu (2009), Paul (2007), Nadler (2008)].

The principal directions of variation correspond to the eigenvectors of the covariance matrix, and in high-dimensions consistent estimation of the eigenvectors is generally not possible without additional assumptions about the covariance matrix or its eigenstructure. Much of the recent development in PCA has focused on methodology that applies the concept of sparsity to the estimation of individual eigenvectors [examples include Jolliffe, Trendafilov and Uddin (2003), d’Aspremont et al. (2007), Zou, Hastie and Tibshirani (2006), Shen and Huang (2008), Witten, Tibshirani and Hastie (2009), Journée et al. (2010)]. Theoretical developments on sparsity and PCA include consistency [Johnstone and Lu (2009), Shen, Shen and Marron (2013)], variable selection properties [Amini and Wainwright (2009)], rates of convergence and minimaxity [Vu and Lei (2012a)], but have primarily been limited to results about estimation of the leading eigenvector. Very recently, Birnbaum et al. (2013) established minimax lower bounds for the estimation of individual eigenvectors. However, an open problem that has remained is whether sparse PCA methods can optimally estimate the subspace spanned by the leading eigenvectors, that is, the principal subspace of variation.

The subspace estimation problem is directly connected to dimension reduction and is important when there may be more than one principal component of interest. Indeed, typical applications of PCA use the projection onto the principal subspace to facilitate exploration and inference of important features of the data. In that case, the assumption that there are distinct principal directions of variation is mathematically convenient but unnatural: it avoids the problem of unidentifiability of eigenvectors by imposing an artifactual choice of principal axes. Dimension reduction by PCA should emphasize subspaces rather than eigenvectors.

In this paper, we study sparse principal subspace estimation in high-dimensions. We present nonasymptotic minimax lower and upper bounds for estimation of both row sparse and column sparse principal subspaces. Our upper bounds are constructive and apply to a wide class of distributions and covariance matrices. In the row sparse case they are optimal up to constant factors, while in the column sparse case they are nearly optimal. As an illustration, one consequence of our results is that the order of the minimax mean squared estimation error of a row sparse $d$-dimensional principal subspace (for $d\ll p$) is

To our knowledge, the only other work that has considered sparse principal subspace estimation is that of Ma (2013). He proposed a sparse principal subspace estimator based on iterative thresholding, and derived its rate of convergence under a spiked covariance model (where the covariance matrix is assumed to be a rank-$d$ perturbation of the identity) similar to that in Birnbaum et al. (2013). He showed that it nearly achieves the optimal rate when estimating a single eigenvector, but was not able to track its dependence on the dimension of the principal subspace.

We obtain the minimax upper bounds by analyzing a sparsity constrained principal subspace estimator and showing that it attains the optimal error (up to constant factors). In comparison to most existing works in the literature, we show that the upper bounds hold without assuming a spiked covariance model. This spiked covariance assumption seems to be necessary for two reasons. The first is that it simplifies analyses and enables the exploitation of special properties of the multivariate Gaussian distribution. The second is that it excludes the possibility of the variables having equal variances. Estimators proposed by Paul (2007), Johnstone and Lu (2009), and Ma (2013) require an initial estimate based on diagonal thresholding—screening out variables with small sample variances. Such an initial estimate will not work when the variables have equal variances or have been standardized. The spiked covariance model excludes that case and, in particular, does not allow PCA on correlation matrices.

A key technical ingredient in our analysis of the subspace estimator is a novel variational form of the Davis–Kahan $\sin\Theta$ theorem (see Corollary 4.1) that may be useful in other regularized spectral estimation problems. It allows us to bound the estimation error using some recent advanced results in empirical process theory, without Gaussian or spiked covariance assumptions. The minimax lower bounds follow the standard Fano method framework [e.g., Yu (1997)], but their proofs involve nontrivial constructions of packing sets in the Stiefel manifold. We develop a generic technique that allows us to convert global packing sets without orthogonality constraints into local packing sets in the Stiefel manifold, followed by a careful combinatorial analysis on the cardinality of the resulting matrix class.

The remainder of the paper is organized as follows. In the next section, we introduce the sparse principal subspace estimation problem and formally describe our minimax framework and estimator. In Section 3, we present our main conditions and results, and provide a brief discussion about their consequences and intuition. Section 4 outlines the key ideas and main steps of the proof. Section 5 concludes the paper with discussion of related problems and practical concerns. Appendices A, B contain the details in proving the lower and upper bounds. The major steps in the proofs require some auxiliary lemmas whose proofs we defer to Appendices C, D.

Subspace estimation

and the orthogonal projector of $\mathcal{S}$ is given by $\Pi_{\mathcal{S}}=VV^{T}$, where $V$ is the $p\times d$ matrix with columns $v_{1},\ldots,v_{d}$.

In practice, $\Sigma$ is unknown, so $\mathcal{S}$ must be estimated from the data. Standard PCA replaces (2) with an empirical version. This leads to the spectral decomposition of the sample covariance matrix

where $\bar{X}$ is the sample mean, and estimating $\mathcal{S}$ by the span of the leading $d$ eigenvectors of $S_{n}$. In high-dimensions however, the eigenvectors of $S_{n}$ can be inconsistent estimators of the eigenvectors of $\Sigma$. Additional structural constraints are necessary for consistent estimation of $\mathcal{S}$.

where $a_{*j}$ denotes the $j$th column of $A$. This is not coordinate-independent. We define the column sparse subspaces to be those that have some orthonormal basis with small $(*,q)$-norm. {definition*}[(Column sparse subspaces)] For $0\leq q<2$ and $1\leq R_{q}\leq p^{1-q/2}$,

2 Parameter space

for $i=1,\ldots,n$, where $\|\cdot\|_{\psi_{\alpha}}$ is the Orlicz $\psi_{\alpha}$-norm [e.g., van der Vaart and Wellner (1996), Chapter 2] defined for $\alpha\geq 1$ as

This ensures that all one-dimensional marginals of $X_{i}$ have sub-Gaussian tails. We also assume that the eigengap $\lambda_{d}-\lambda_{d+1}>0$ so that the principal subspace $\mathcal{S}$ is well defined. Intuitively, $\mathcal{S}$ is harder to estimate when the eigengap is small. This is made precise by the effective noise variance

It turns out that this is a key quantity in the estimation of $\mathcal{S}$, and that it is analogous to the noise variance in linear regression. Let

denote the class of distributions on $X_{1},\ldots,X_{n}$ that satisfy (4), $\sigma_{d}^{2}\leq\sigma^{2}$, and $\mathcal{S}\in\mathcal{M}_{q}(R_{q})$. Similarly, let

denote the class of distributions that satisfy (4), $\sigma_{d}^{2}\leq\sigma^{2}$, and $\mathcal{S}\in\mathcal{M}_{q}^{*}(R_{q})$.

3 Subspace distance

for $k=1,\ldots,d$ and the angle operator between $\mathcal{E}$ and $\mathcal{F}$ is the $d\times d$ matrix

In other words, $EF^{\perp}$ has at most $d$ nonzero singular values and the nonzero singular values of $E-F$ are the nonzero singular values of $EF^{\perp}$, each counted twice.

We will frequently use these identities. For simplicity, we will overload notation and write

In other words, the distance between two subspaces is equivalent to the minimal distance between their orthonormal bases.

4 Sparse subspace estimators

Here we introduce an estimator that achieves the optimal (up to a constant factor) minimax error for row sparse subspace estimation. To estimate a row sparse subspace, it is natural to consider the empirical minimization problem corresponding to (2) with an additional sparsity constraint corresponding to $\mathcal{M}_{q}(R_{q})$.

We define the row sparse principal subspace estimator to be a solution of the following constrained optimization problem:

For our analysis, it is more convenient to work on the Stiefel manifold. Let $\langle A,B\rangle:=\operatorname{trace}(A^{T}B)$ for matrices $A,B$ of compatible dimension. It is straightforward to show that following optimization problem is equivalent to (2.4):

If $\hat{V}$ is a global maximizer of (2.4), then $\operatorname{col}(\hat{V})$ is a solution of (2.4). When $q=1$, the estimator defined by (2.4) is essentially a generalization to subspaces of the Lasso-type sparse PCA estimator proposed by Jolliffe, Trendafilov and Uddin (2003). A similar idea has also been used by Chen, Zou and Cook (2010) in the context of sufficient dimension reduction. The constraint set in (2.4) is clearly nonconvex, however this is unimportant, because the objective function is convex and we know that the maximum of a convex function over a set $D$ is unaltered if we replace $D$ by its convex hull. Thus, (2.4) is equivalent to a convex maximization problem. Finding a global maximum of convex maximization problems is computationally challenging and efficient algorithms remain to be developed. Nevertheless, in the most popular case $q=1$, some algorithms have been proposed with promising empirical performance [Shen and Huang (2008), Witten, Tibshirani and Hastie (2009)].

We define the column sparse principal subspace estimator analogously to the row sparse principal subspace estimator, using the column sparse subspaces $\mathcal{M}_{q}^{*}(R_{q})$ instead of the row sparse ones. This leads to the following equivalent Grassmann and Stiefel manifold optimization problems:

Main results

In this section, we present our main results on the minimax lower and upper bounds on sparse principal subspace estimation over the row sparse and column sparse classes.

To highlight the key results with minimal assumptions, we will first consider the simplest case where $q=0$. Consider the following two conditions.

$4\leq p-d$ and $2d\leq R_{q}-d\leq(p-d)^{1-{q/2}}$.

Condition 1 is necessary for the existence of a consistent estimator (see Theorems A.1 and A.2). Without Condition 1, the statements of our results would be complicated by multiple cases to deal with the fact that the subspace distance is bounded above by $\sqrt{d}$. The lower bounds on $p-d$ and $R_{q}-d$ are minor technical conditions that ensure our nonasymptotic bounds are nontrivial. Similarly, the upper bound on $R_{q}-d$ is only violated in trivial cases (detailed discussion given below).

Here, as well as in the entire paper, $c$ denotes a universal, positive constant, not necessarily the same at each occurrence. This lower bound result reflects two separate aspects of the estimation problem: variable selection and parameter estimation after variable selection. Variable selection refers to finding the variables that generate the principal subspace, while estimation refers to estimating the subspace after selecting the variables. For each variable, we accumulate two types of errors: one proportional to $d$ that reflects the coordinates of the variable in the $d$-dimensional subspace, and one proportional to $\log[(p-d)/(R_{0}-d)]$ that reflects the cost of searching for the $R_{0}$ active variables. We prove Theorem 3.1 in Appendix A.

The interpretation for these two quantities is natural. First, $T$ measures the relative sparsity of the problem. Roughly speaking, it ranges between and $1$ when the sparsity constraint in (2.4) is active, though the “sparse” regime generally corresponds to $T\ll 1$. The second quantity, $\gamma$ corresponds to the classic mean squared error (MSE) of standard PCA. The problem is low-dimensional if $\gamma$ is small compared to $T$. We impose the following condition to preclude this case.

There is a constant $a<1$ such that $T^{a}\leq\gamma^{q/2}$.

This condition lower bounds the classic MSE in terms of the sparsity and is mild in high-dimensional situations. When $a=q/2$, for example, Condition 3 reduces to

Let $q\in(0,2)$. If Conditions 1 to 3 hold, then

This result generalizes Theorem 3.1 and reflects the same combination of variable selection and parameter estimation. When Condition 3 does not hold, the problem is outside of the sparse, high-dimensional regime. As we show in the proof, there is actually a “phase transition regime” between the high-dimensional sparse and the classic dense regimes for which sharp minimax rate remains unknown. A similar phenomenon has been observed in Birnbaum et al. (2013).

2 Row sparse upper bound

Our upper bound results are obtained by analyzing the estimators given in Section 2.4. The case where $q=0$ is the clearest, and we begin by stating a weaker, but simpler minimax upper bound for the row sparse class.

Let $\hat{\mathcal{S}}$ be any global maximizer of (2.4). If $6\sqrt{R_{0}(d+\log p)}\leq\sqrt{n}$, then

Although (2.4) may not have a unique global optimum, Theorem 3.3 shows that any global optimum will be within a certain radius of the principal subspace $\mathcal{S}$. The proof of Theorem 3.3, given in Section 4.2, is relatively simple but still nontrivial. It also serves as a prototype for the much more involved proof of our main upper bound result stated in Theorem 3.4 below. We note that the rate given by Theorem 3.3 is off by a $\lambda_{1}/\lambda_{d+1}$ factor that is due to the specific approach taken to control an empirical process in our proof of Theorem 3.3.

To state the main upper bound result with optimal dependence on ($n$, $p$, $d$, $R_{q}$, $\sigma^{2}$), we first describe some regularity conditions. Let

where $c_{1}$ and $c_{3}$ are positive constants involved in the empirical process arguments. Equations (12) to (15) require that $\varepsilon_{n}$, the minimax rate of estimation (except the factor involving $\lambda$), to be small enough, compared to empirical process constants and some polynomials of $\lambda$. Such conditions are mild in the high dimensional, sparse regime, since to some extent, they are qualitatively similar and analogous to Conditions 1 to 3 required by the lower bound.

Conditions (12) to (15) are general enough to allow $R_{q}$, $d$ and $\lambda_{j}$ ($j=1,d,d+1$) to scale with $n$. For example, consider the case $q=0$, and let $d=n^{a}$, $R_{0}=n^{b}$, $p=n^{c}$, $\lambda_{1}=n^{r_{1}}$, $\lambda_{d}=n^{r_{2}}$, $\lambda_{d+1}=n^{r_{3}}$, where $0<a<b<c$, and $r_{1}\geq r_{2}>r_{3}$. Note that the $r_{j}$’s can be negative. Then it is straightforward to verify that conditions (12) to (15) hold for large values of $n$ whenever $a+b<1$ and $r_{1}<r_{2}+(1-a)/2$. Condition (12) implies that $d$ cannot grow faster than $\sqrt{n}$.

with probability at least $1-4/(n-1)-6\log n/n-p^{-1}$.

Theorem 3.4 is presented in terms of a probability bound instead of an expectation bound. This stems from technical aspects of our proof that involve bounding the supremum of an empirical process over a set of random diameter. For $q\in$, the upper bound matches our lower bounds (Theorems 3.1 and 3.2) for the entire tuple ($n$, $p$, $d$, $R_{q}$, $\sigma^{2}$) up to a constant if

for some constant $c<1$. To see this, combining this additional condition and Condition 2, the term $\log\frac{p^{1-q/2}}{R_{q}}$ in the lower bound given in Theorem 3.2 is within a constant factor of $\log p$ in the upper bound given in Theorem 3.4. It is straightforward to check that the other terms in lower and upper bounds agree up to constants with obvious correspondence. Moreover, we note that the additional condition (16) is only slightly stronger than the last inequality in Condition 2. The proof of Theorem 3.4 is in Appendix B.1.

Using the probability upper bound result and the fact that $\|{\sin\Theta}(\hat{\mathcal{S}},\mathcal{S})\|_{F}^{2}\leq d$, one can derive an upper bound in expectation.

Under the same condition as in Theorem 3.4, we have for some constant $c$,

The expectation upper bound has an additional $d(\log n/n+1/p)$ term that can be further reduced by refining the argument (see Remark 3 below). It is not obvious if one can completely avoid such a term. But in many situations it is dominated by the first term. Again, we invoke the scaling considered in Remark 1. When $q=0$, the first term is of order $n^{a+b+(r_{1}+r_{3})/2-r_{2}-1}$, and the additional term is $n^{a-1}\log n+n^{a-c}$, which is asymptotically negligible if $b>r_{2}-(r_{1}+r_{3})/2+(1-c)_{+}$.

Given any $r>0$, it is easy to modify the proof of Theorem 3.4 [as well as conditions (12) to (15)] such that the results of Theorem 3.4 and Corollary 3.1 hold with $c$ replaced by some constant $c(r)$, and the probability bound becomes $1-4/(n^{r}-1)-6\log n/n^{r}-1/p^{r}$.

3 Column sparse lower bound

By modifying the proofs of Theorems 3.1 and 3.2, we can obtain lower bound results for the column sparse case that are parallel to the row sparse case. For brevity, we present the $q=0$ and $q>0$ cases together. The analog of $T$, the degree of sparsity, for the column sparse case is

and the analogs of Conditions 2 and 3 are the following.

$4d\leq p-d$ and $d\leq d(R_{q}-1)\leq(p-d)^{1-{q/2}}$.

There is a constant $a<1$ such that $T_{*}^{a}\leq\gamma^{q/2}$.

Let $q\in[0,2)$. If Conditions 4 and 5 hold, then

4 Column sparse upper bound

A specific challenge in analyzing the column sparse principal subspace problem (2.4) is to bound the supremum of the empirical process

indexed by all $U\in\mathcal{U}(p,d,R_{q},\varepsilon)$ where

Unlike the row sparse matrices, the matrices $UU^{T}$ and $VV^{T}$ are no longer column sparse with the same radius $R_{q}$.

By observing that $\mathcal{M}_{q}^{*}(R_{q})\subseteq\mathcal{M}_{q}(dR_{q})$, we can reuse the proof of Theorem 3.4 to derive the following upper bound for the column sparse class.

with probability at least $1-4/(n-1)-6\log n/n-p^{-1}$.

Corollary 3.2 is slightly weaker than the corresponding result for the row sparse class. It matches the lower bound in Theorem 3.5 up to a constant if

for some constant $c<1$, and $d<C\log p$ for some other constant $C$.

5 A conjecture for the column sparse case

Note that Theorem 3.5 and Corollary 3.2 only match when $d\leq C\log p$. For larger values of $d$, we believe that the lower bound in Theorem 3.5 is optimal and the upper bound can be improved.

[(Minimax error bound for column sparse case)] Under the same conditions as in Corollary 3.2, there exists an estimator $\hat{\mathcal{S}}$ such that

with high probability. As a result, the optimal minimax lower and upper bounds for this case shall be

One reason for the conjecture is based on the following intuition. Suppose that $\lambda_{1}>\lambda_{2}>\cdots>\lambda_{d}>\lambda_{d+1}$ (there is enough gap between the leading eigenvalues) one can recover the individual leading eigenvectors with an error rate whose dependence on $(n,R_{q},p)$ is the same as in the lower bound [cf. Vu and Lei (2012a), Birnbaum et al. (2013)]. As a result, the estimator $\hat{V}=(\hat{v}_{1},\hat{v}_{2},\ldots,\hat{v}_{d})$ shall give the desired upper bound. On the other hand, it remains open to us whether the estimator in (2.4) can achieve this rate for $d$ much larger than $\log p$.

Sketch of proofs

For simplicity, we focus on the row sparse case with $q=0$, assuming also the high dimensional and sparse regime. For more general cases, see Theorems A.1 and A.2 in Appendix A.

Our proof of the lower bound features a combination of the general framework of the Fano method and a careful combinatorial analysis of packing sets of various classes of sparse matrices. The particular challenge is to construct a rich packing set of the parameter space $\mathcal{P}_{q}(\sigma^{2},R_{q})$. We will consider centered $p$-dimensional Gaussian distributions with covariance matrix $\Sigma$ given by

for $0\leq\varepsilon\leq 1$. We have the following generic method for lower bounding the minimax risk of estimating the principal subspace of a covariance matrix. It is proved in Appendix A as a consequence of Lemmas A.1 to A.3.

then every estimator $\hat{\mathcal{A}}$ of $\mathcal{A}_{i}:=\operatorname{col}(A_{\varepsilon}(J_{i}))$ satisfies

Note that if $\|J\|_{2,0}\leq R_{0}-d$, then $\|A_{\varepsilon}(J)\|_{2,0}\leq R$. Thus Lemma 4.1 with appropriate choices of $J_{i}$ can yield minimax lower bounds over $p$-dimensional Gaussian distributions whose principal subspace is $R_{0}$ row sparse.

This yields a minimax lower bound that reflects the complexity of post-selection estimation. Putting these two results together, we have for a subset of Gaussian distributions $G\subseteq\mathcal{P}_{0}(\sigma^{2},R_{q})$ the minimax lower bound:

2 The upper bound

The upper bound proof requires a careful analysis of the behavior of the empirical maximizer of the PCA problem under sparsity constraints. The first key ingredient is to provide a lower bound of the curvature of the objective function at its global maxima. Traditional results of this kind, such as Davis–Kahan sin$\Theta$ theorem and Weyl’s inequality, are not sufficient for our purpose.

The following lemma, despite its elementary form, has not been seen in the literature (to our knowledge). It gives us the right tool to bound the curvature of the matrix functional $F\mapsto\langle A,F\rangle$ at its point of maximum on the Grassmann manifold.

Lemma 4.2 is proved in Appendix C.2. An immediate corollary is the following alternative to the traditional matrix perturbation approach to bounding subspace distances using the Davis–Kahan $\sin\Theta$ theorem and Weyl’s inequality.

In addition to the hypotheses of Lemma 4.2, if $B$ is a symmetric matrix and $F$ satisfies

The corollary is different from the Davis–Kahan $\sin\Theta$ theorem because the orthogonal projector $F$ does not have to correspond to a subspace spanned by eigenvectors of $B$. $F$ only has to satisfy

This condition is suited ideally for analyzing solutions of regularized and/or constrained maximization problems where $E$ and $F$ are feasible, but $F$ is optimal. In the simplest case, where $g\equiv 0$, combining (21) with the Cauchy–Schwarz inequality and (6) recovers a form of the Davis–Kahan $\sin\Theta$ theorem in the Frobenius norm:

Applying Corollary 4.1 with $B=S_{n}$, $A=\Sigma$, $E=VV^{T}$, $F=\hat{V}\hat{V}^{T}$, and $g\equiv 0$, we have

Obtaining a sharp upper bound for $\langle S-\Sigma,\hat{V}\hat{V}^{T}-VV^{T}\rangle$ is nontrivial. First, one needs to control $\sup_{F\in\mathcal{F}}\langle S-\Sigma,F\rangle$ for some class $\mathcal{F}$ of sparse and symmetric matrices. This requires some results on quadratic form empirical process. Second, in order to obtain better bounds, we need to take advantage of the fact that $\hat{V}\hat{V}^{T}-VV^{T}$ is probably small. Thus, we need to use a peeling argument to deal with the case where $\mathcal{F}$ has a random (but probably) small diameter. These details are given in Appendices B.1 and D. Here we present a short proof of Theorem 3.3 to illustrate the idea.

because $\|\hat{V}\hat{V}^{T}-VV^{T}\|_{F}^{2}=2\hat{\varepsilon}^{2}$ by (6). Let

The empirical process $\langle S_{n}-\Sigma,UU^{T}\rangle$ indexed by $U$ is a generalized quadratic form, and a sharp bound of its supremum involves some recent advances in empirical process theory due to Mendelson (2010) and extensions of his results. By Corollary 4.1 of Vu and Lei (2012b), we have

because $U\in\mathcal{U}(R_{0})$. Using a standard $\delta$-net argument (see Propositions D.1 and D.2), we have, when $p>5$,

The proof is complete since we assume that $6\sqrt{R_{0}(d+\log p)}\leq\sqrt{n}$.

Discussion

There is a natural correspondence between the sparse principal subspace optimization problem (2.4) and some optimization problems considered in the sparse regression literature. We have also found that there is a correspondence between minimax results for sparse regression and those that we presented in this article. In spite of these connections, results on computation for sparse principal subspaces (and sparse PCA) are far less developed than for sparse regression. In this final section, we will discuss the connections with sparse regression, both optimization and minimax theory, and then conclude with some open problems for sparse principal subspaces.

for $0<q\leq 1$ and similarly for $q=0$. When $q=1$, this corresponds to a “group Lasso” penalty where entries in the same row of $U$ are penalized simultaneously [Yuan and Lin (2006), Zhao, Rocha and Yu (2009)]. The idea being that as $\tau_{q}$ varies, a variable should enter/exit all $d$ coordinates simultaneously. In the column sparse case, when $q=1$ the analogous penalized multivariate regression problem has a penalty which encourages each column of $U$ to be sparse, but does not require that the pattern of sparsity to be the same across columns.

agreeing with our minimax lower and upper bounds for the row sparse principal subspace problem.

2 Practical concerns

Although the minimax optimal estimators that we propose do not require knowledge of the noise-to-signal ratio $\sigma^{2}$, they do require knowledge of (or an upper bound on) the sparsity $R_{q}$. It is not hard to modify our techniques to produce an estimator that gives up adaptivity to $\sigma^{2}$ in exchange for adaptivity to $R_{q}$. One could do this by using penalized versions of our estimators with a penalty factor proportional to $\sigma^{2}$. An extension along this line has already been considered by Lounici (2013) for the $d=1$ case. A more interesting question is whether or not there exist fully adaptive principal subspace estimators.

Under what conditions can one find an estimator that achieves the minimax optimal error without requiring knowledge of either $\sigma^{2}$ or $R_{q}$? Works by Paul (2007) and Ma (2013) on refinements of diagonal thresholding for the spiked covariance model seems promising on this front, but as we mentioned in the Introduction, the spiked covariance model is restrictive and necessarily excludes the common practice of standardizing variables. Is it possible to be adaptive outside the spiked covariance model? One possible approach can be described in the following three steps. (1) use a conservative choice of $R_{q}$ (say, $p^{a}$, for some $0<a<1$); (2) estimate $\sigma^{2}$ using eigenvalues obtained from the sparsity constrained principal subspace estimator; and (3) use a sparsity penalized principal subspace estimator with $\sigma^{2}$ replaced by its estimate. We will pursue this idea in further detail in future work.

Appendix A Lower bound proofs

Theorems 3.1, 3.2 and 3.5 are consequences of three more general results stated below. An essential part of the strategy of our proof is to analyze the variable selection and estimation aspects of the problem separately. We will consider two types of subsets of the parameter space that capture the essential difficulty of each aspect: one where the subspaces vary over different subsets of variables, and another where the subspaces vary over a fixed subset of variables. The first two results give lower bounds for each aspect in the row sparse case. Theorems 3.1 and 3.2 follow easily from them. The third result directly addresses the proof of Theorem 3.5.

Let $q\in[0,2)$ and $(p,d,R_{q})$ satisfy

There exists a universal constant $c>0$ such that every estimator $\hat{\mathcal{S}}$ satisfies the following. If $T<\gamma^{q/2}$, then

The case $q=0$ is particularly simple, because $T<\gamma^{q/2}=1$ holds trivially. In that case, Theorem A.1 asserts that

When $q\in(0,2)$ the transition between the $T<\gamma^{q/2}$ and $T\geq\gamma^{q/2}$ regimes involves lower order ($\log\log$) terms that can be seen in (A.2). Under Condition 3, (A.1) can be simplified to

Let $q\in[0,2)$ and $(p,d,R_{q})$ satisfy

and let $T$ and $\gamma$ be defined as in (11). There exists a universal constant $c>0$ such that every estimator $\hat{\mathcal{S}}$ satisfies the following. If $T<(d\gamma)^{q/2}$, then

This result with (A) implies Theorem 3.1, and with (A) it implies Theorem 3.2.

Let $q\in[0,2)$ and $(p,d,R_{q})$ satisfy

and recall the definition of $T_{*}$ in (17). There exists a universal constant $c>0$ such that every estimator $\hat{\mathcal{S}}$ satisfies the following. If $T_{*}<\gamma^{q/2}$, then

In the next section we setup a general technique, using Fano’s inequality and Stiefel manifold embeddings, for obtaining minimax lower bounds in principal subspace estimation problems. Then we move on to proving Theorems A.1 and A.3.

Our main tool for proving minimax lower bounds is the generalized Fano method. We quote the following version from Yu (1997), Lemma 3.

and, the Kullback–Leibler (KL) divergence

Then every $\mathcal{A}$-measurable estimator $\hat{\theta}$ satisfies

where $b>0$. The noise-to-signal ratio of the principal $d$-dimensional subspace of these covariance matrices is

and can choose $b$ to achieve any $\sigma^{2}>0$. The KL divergence between these multivariate Normal distributions has a simple, exact expression given in the following lemma. The proof is straightforward and contained in Appendix C.1.

By using Lemma A.3 in conjunction with Lemmas A.1 and A.2, we have the following generic method for lower bounding the minimax risk of estimating the principal subspace of a covariance matrix.

then every estimator $\hat{\mathcal{A}}$ of $\mathcal{A}_{i}:=\operatorname{col}(A_{\varepsilon}(J_{i}))$ satisfies

A.2 Proofs of the main lower bounds

Proof of Theorem A.1 The following lemma, derived from Massart [(2007), Lemma 4.10], allows us to analyze the variable selection aspect.

$\|J_{i}-J_{j}\|_{2}^{2}\geq 1/4$ for all $i\neq j$, and

$\log N\geq\max\{cs[1+\log(m/s)],\log(m)\}$, where $c>1/30$ is an absolute constant.

Applying Lemma A.4, with $k=1$, $\delta_{N}=1/2$, and $b$ chosen so that $(1+b)/b^{2}=\sigma^{2}$, yields

If we can choose $\rho\in(0,1]$ such that (37) is satisfied, then by (A.2),

Choose $\rho\in(0,1]$ to be the unique solution of the equation

We will verify that $\varepsilon$ and $\rho$ satisfy (37). The assumption that $1\leq R_{q}-d$ guarantees that $\varepsilon^{2q}\leq(R_{q}-d)^{2}$, because $\varepsilon^{2q}\leq 1$. If $T<\gamma^{q/2}$, then

If $T\geq\gamma^{q/2}$, then $\rho=1$ and

Now we substitute (38) and the definitions of $\gamma$ and $T$ into the above inequality to get the following lower bounds. If $T<\gamma^{q/2}$, then

If $T\geq\gamma^{q/2}$, then $\gamma\rho(1-\log\rho)=\gamma$ and

$\|\sin(J_{i},J_{j})\|_{F}\geq\sqrt{k}\delta$ for all $i\neq j$, and

$\log N\geq k(s-k)\log(c_{2}/\delta)$, where $c_{2}>0$ is an absolute constant.

To apply this result to Lemma A.4 we will use Proposition 2.2 to convert the lower bound on the subspace distance into a lower bound on the Frobenius distance between orthonormal matrices. Thus,

for all $i\neq j$. The rest of this proof mirrors that of Theorem A.1. Let $\varepsilon\in[0,1/\sqrt{2}]$ and apply Lemma A.4 to get

So $\varepsilon$ and $\rho$ must satisfy the constraint

where the right-hand side is an assumption of the lemma. That verifies one of the inequalities in (42). If $T<(d\gamma)^{q/2}$, then

If $T\geq(d\gamma)^{q/2}$, then $\rho=1$ and

Finally, we substitute the definition of $\gamma$ and (44) into the above inequality to get the following lower bounds. If $T<(d\gamma)^{q/2}$, then

Proof of Theorem A.3 The proof is a modification of the proof of Theorem A.1. The difficulty of the problem is captured by the difficulty of variable selection within each column of $V$. Instead of using a single hypercube construction as in the proof of Theorem A.1, we apply a hypercube construction on each of the $d$ columns. We do this by dividing the $(p-d)\times d$ matrix into $d$ submatrices of size $\lfloor(p-d)/d\rfloor\times d$, that is, constructing matrices of the form

and confining the hypercube construction to the $k$th column of each $\lfloor(p-d)/d\rfloor\times d$ matrix $B_{k}$, $k=1,\ldots,d$. This ensures that the resulting $(p-d)\times d$ matrix has orthonormal columns with disjoint supports.

$\|J_{i}-J_{j}\|_{2}^{2}\geq 1/4$ for all $i\neq j$, and

$\log M\geq\max\{cs(1+\log(m/s)),\log m\}$, where $c>1/30$ is an absolute constant.

$\|H\|_{*,0}\leq s$ for all $H\in\mathcal{H}^{s}$.

$\|H_{1}-H_{2}\|_{2}^{2}\geq d/8$ for all $H_{1},H_{2}\in\mathcal{H}^{s}$ such that $H_{1}\neq H_{2}$.

$\log N:=\log|\mathcal{H}^{s}|\geq\max\{cds(1+\log(m/s)),\log m\}$, where $c>0$ is an absolute constant.

Note that the lower bound of $\log m$ in the third item arises by considering the packing set whose $N$ elements consist of matrices whose columns in $B_{1},\ldots,B_{d}$ are all equal to some $J_{i}$ for $i=1,\ldots,M$. This ensures that $\log N\geq\log M\geq\log m$. From here, the proof is a straightforward modification of proof of Theorem A.1 with the substitution of $p-d$ by $(p-d)/d$. For brevity we will only outline the major steps.

Recall the definitions of $T_{*}$ and $\gamma$ in (17). Apply Lemma A.4 with the subset $\mathcal{H}^{s}$, $k=d$, $\delta_{N}=\sqrt{d}/\sqrt{8}$, and $b$ chosen so that $(1+b)/b^{2}=\sigma^{2}$. Then

by the assumption that $(p-d)/d\geq 4$, and

where $c_{1}>0$ is a sufficiently small constant, the assumption that $d<\break d(R_{q}-1)$, and letting $\rho$ be the unique solution of the equation

We conclude that every estimator $\hat{V}$ satisfies

and we have the following explicit lower bounds. If $T_{*}<\gamma^{q/2}$, then

Appendix B Upper bound proofs

$\Sigma$ and $S_{n}$ are both invariant under translations of $\mu$. Since our estimators only depend on $X_{1},\ldots,X_{n}$ only through $S_{n}$, we will assume without loss of generality that $\mu=0$ for the remainder of the paper. The sample covariance matrix can be written as

It can be show that $\bar{X}\bar{X}^{T}$ is a higher order term that is negligible [see the proofs in Vu and Lei (2012a), for an example of such arguments]. Therefore, we will ignore this term and focus on the dominating $\frac{1}{n}\sum_{i=1}^{n}X_{i}X_{i}^{T}$ term in our proofs below. {pf*}Proof of Theorem 3.4 Again, we start from Corollary 4.1, which gives

To get the correct dependence on $\lambda_{i}$ and for general values of $q$, we need a more refined analysis to control the random variable $\langle S_{n}-\Sigma,\hat{V}\hat{V}^{T}-VV^{T}\rangle$. Let

Recall that for an orthogonal projector $\Pi$ we write $\Pi^{\perp}:=I-\Pi$. By Proposition C.1 we have

We will control $T_{1}$ (the upper-quadratic term), $T_{2}$ (the cross-product term), and $T_{3}$ (the lower-quadratic term) separately.

where $c_{1}$ is a universal constant. Define

Combining (B.1) and (B.1) we obtain, for all $t>0$, $0<q<1$,

The case where $q=0$ is simpler and omitted. Now define

Taking $t=t_{2,2}$ in (52) and using the tail bound result in Lemma D.1, we have

The bound on $T_{3}$ involves a quadratic form empirical process over a random set. Let $\varepsilon\geq 0$ and define

Then by Lemma D.4, we have, with some universal constants $c_{3}$, for $x>0$

Let $T_{3}(U)=\langle W,\Pi^{\perp}UU^{T}\Pi^{\perp}\rangle$, for all $U\in\mathcal{U}_{p}(R_{q})$, where

Define function $g(\varepsilon)=\varepsilon_{n}\varepsilon^{2}+\varepsilon_{n}^{2}\varepsilon+\varepsilon_{n}^{4}$. Then for all $\varepsilon\geq 0$, we have $g(\varepsilon)\geq\varepsilon_{n}^{4}\geq 4d^{3}/n^{2}$. On the other hand, if $\varepsilon=\|{\sin\Theta}(U,V)\|_{F}$, then $\varepsilon^{2}\leq 2d$ and hence $g(\varepsilon)\leq g(\sqrt{2d})=2d+\sqrt{2d}+1$. Let $\mu=\varepsilon_{n}^{4}$ and $J=\lceil\log_{2}(g(\sqrt{2d})/\mu)\rceil$. Then we have $J\leq 3\log n+6/5$.

Note that $g$ is strictly increasing on $[0,\sqrt{2d}]$. Then we have the following peeling argument:

Putting things together

Now recall the conditions in (12) to (15). On $\Omega_{1}^{c}\cap\Omega_{2}^{c}\cap\Omega_{3}^{c}$, we have, from (45) that

Appendix C Additional proofs

Proof of Proposition 2.2 Let $\gamma_{i}$ be the cosine of the $i$th canonical angle between the subspaces spanned by $V_{1}$ and $V_{2}$. By Theorem II.4.11 of Stewart and Sun (1990),

Apply the trigonometric identity $\sin^{2}\theta=1-\cos^{2}\theta$ to the preceding display to conclude the proof.

Proof of Lemma A.2 Write $\Sigma_{i}=\Sigma(A_{i})$ for $i=1,2$. Since $\Sigma_{1}$ and $\Sigma_{2}$ are nonsingular and have the same determinant,

Proof of Lemma A.3 By Proposition 2.1 and the definition of $A_{\varepsilon}(\cdot)$,

The upper bound follows from Proposition 2.2:

Proof of Lemma A.5 Let $s_{0}=\lfloor\min(m/e,s)\rfloor$. The assumptions that $m/e\geq 1$ and $s\geq 1$ guarantee that $s_{0}\geq 1$. According to Massart [(2007), Lemma 4.10] [with $\alpha=7/8$ and $\beta=8/(7e)$], there exists a subset $\Omega_{m}^{s_{0}}\subseteq\{0,1\}^{m}$ satisfying the following properties:

$\|\omega\|_{0}=s_{0}$ for all $\omega\in\Omega_{m}^{s_{0}}$,

$\|\omega-\omega^{\prime}\|_{0}>s_{0}/4$ for all distinct pairs $\omega,\omega^{\prime}\in\Omega_{m}^{s_{0}}$, and

$\log|\Omega_{m}^{s_{0}}|\geq cs_{0}\log(m/s_{0})$, where $c>0.251$.

The cardinality of $\{J_{1},\ldots,J_{N}\}$ satisfies

As a function of $s_{0}$, the above right-hand side is increasing on the interval $[0,m/e]$. Since $\min(m/e,s)/2\leq s_{0}$ belongs to that interval

where $(c/2)(1+e)^{-1}>1/30$. If the above right-hand side is $\leq\log m$, then we may repeat the entire argument from the beginning with $\{J_{1},\ldots,J_{N}\}$ taken to be the $N=m$ vectors $\{(1,0,\ldots,0),(0,1,0,\ldots,0),\ldots,(0,\ldots,0,1)\}\subseteq\{0,1\}^{m}$. That yields, in combination with (54),

C.2 Proofs related to the upper bounds

Proof of Lemma 4.2 For brevity, denote the eigenvalues of $A$ by $\lambda_{d}:=\lambda_{d}(A)$. Let $A=\sum_{i=1}^{p}\lambda_{i}u_{i}u_{i}^{T}$ be the spectral decomposition of $A$ so that $E=\sum_{i=1}^{d}u_{i}u_{i}^{T}$ and $E^{\perp}=\sum_{i=d+1}^{p}u_{i}u_{i}^{T}$. Then

Since orthogonal projectors are idempotent,

Now apply Proposition 2.1 to conclude that

If $W$ is symmetric, and $E$ and $F$ are orthogonal projectors, then

and the symmetry of $W$, $F$ and $E$, we can write

Appendix D Empirical process related proofs

This section is dedicated to proving the following bound on the cross-product term.

There exists a universal constant $c>0$ such that

The proof of Lemma D.1 builds on the following two lemmas. They are adapted from Lemmas 2.2.10 and 2.2.11 of van der Vaart and Wellner (1996).

Let $Y_{1},\ldots,Y_{n}$ be independent random variables with zero mean. Then

Let $Y_{1},\ldots,Y_{m}$ be arbitrary random variables that satisfy the bound

for all $t>0$ (and $i$) and fixed $a,b>0$. Then

We bound $\|\Pi^{\perp}(S_{n}-\Sigma)\Pi\|_{2,\infty}$ by a standard $\delta$-net argument.

and $\|u_{*}-u\|_{2}\leq\delta$. Then by the Cauchy–Schwarz inequality,

The following bound on the covering number of the sphere is well known [see, e.g., Ledoux (2001), Lemma 3.18].

Let $X$ and $Y$ be random variables. Then

Let $A=X/\|X\|_{\psi_{2}}$ and $Y/\|Y\|_{\psi_{2}}$. Using the elementary inequality

Multiplying both sides of the inequality by $\|X\|_{\psi_{2}}\|Y\|_{\psi_{2}}$ gives the desired result.

where $e_{j}$ is the $j$th column of $I_{p\times p}$. Taking $\delta=1/2$, by Proposition D.2 we have $|N_{\delta}|\leq 5^{d}$.

is the sum of independent random variables with mean zero. By Proposition D.3, the summands satisfy

Recall that $\|Z\|_{\psi_{2}}^{2}=1$. Then Bernstein’s inequality (Lemma D.2) implies that for all $t>0$ and every $u\in\mathcal{N}_{\delta}$

D.2 The quadratic terms

Let $\varepsilon\geq 0$, $q\in(0,1]$, and

There exist constants $c>0$ and $c_{1}$ such that for all $x\geq c_{1}$,

and $\mathcal{Z}$ is a $p\times d$ matrix with i.i.d. $\mathcal{N}(0,1)$ entries. As a consequence, we have, for another constant $c_{2}$

Moreover, we have, for another numerical constant $c^{\prime}$,

The first part follows from Corollary 4.1 of Vu and Lei (2012b). It remains for us to prove the “moreover” part. By the duality of the $(2,1)$- and $(2,\infty)$-norms,

By (24) and the fact that the Orlicz $\psi_{2}$-norm bounds the expectation,

Using a similar argument as in the proof of Lemma D.1, for all $t>0$ and every $u\in\mathcal{N}_{\delta}$

where $\sigma=2\|Z_{1}\|_{\psi_{2}}^{2}\lambda_{1}$. Then Lemma D.3 implies that

where $C>0$ is a constant. Choosing $\delta=1/3$ and applying Proposition D.2 yields $|\mathcal{N}_{\delta}|\leq 7^{d}$ and

Acknowledgments

The research reported in this article was completed while V. Q. Vu was visiting the Department of Statistics at Carnegie Mellon University. He thanks them for their hospitality and support. We also thank the referees for their helpful comments.

References