Minimax Rates of Estimation for Sparse PCA in High Dimensions

Introduction

High-dimensional data problems, where the number of variables $p$ exceeds the number of observations $n$, are pervasive in modern applications of statistical inference and machine learning. Such problems have increased the necessity of dimensionality reduction for both statistical and computational reasons. In some applications, dimensionality reduction is the end goal, while in others it is just an intermediate step in the analysis stream. In either case, dimensionality reduction is usually data-dependent and so the limited sample size and noise may have an adverse affect. Principal components analysis (PCA) is perhaps one of the most well known and widely used techniques for unsupervised dimensionality reduction. However, in the high-dimensional situation, where $p/n$ does not tend to 0 as $n\to\infty$, PCA may not give consistent estimates of eigenvalues and eigenvectors of the population covariance matrix . To remedy this situation, sparsity constraints on estimates of the leading eigenvectors have been proposed and shown to perform well in various applications. In this paper we prove optimal minimax error bounds for sparse PCA when the leading eigenvector is sparse.

[[, see]Chapter 7.2.3 for example]Izenman:2008. The optimal subspace is determined by spectral decomposition of the population covariance matrix

In practice, $\Sigma$ is not known and so $\Theta$ must be estimated from the data. In that case we replace $\Theta$ by an estimate $\hat{\Theta}$ and reduce the dimension of the data by the mapping $x\mapsto\hat{\Pi}x$, where $\hat{\Pi}=\hat{\Theta}\hat{\Theta}^{T}$. PCA uses the spectral decomposition of the sample covariance matrix

where $\bar{X}$ is the sample mean, and $l_{j}$ and $u_{j}$ are eigenvalues and eigenvectors of $S$ defined analogously to eq. 2. It reduces the dimension of the data to $k$ by the mapping $x\mapsto UU^{T}x$, where $U=(u_{1},\ldots,u_{k})$.

In the classical regime where $p$ is fixed and $n\to\infty$, PCA is a consistent estimator of the population eigenvectors. However, this scaling is not appropriate for modern applications where $p$ is comparable to or larger than $n$. In that case, it has been observed that if $p,n\to\infty$ and $p/n\to c>0$, then PCA can be an inconsistent estimator in the sense that the angle between $u_{1}$ and $\theta_{1}$ can remain bounded away from even as $n\to\infty$.

2 Sparsity Constraints

Estimation in high-dimensions may be beyond hope without additional structural constraints. In addition to making estimation feasible, these structural constraints may also enhance interpretability of the estimators. One important example of this is sparsity. The notion of sparsity is that a few variables have large effects, while most others are negligible. This type of assumption is often reasonable in applications and is now widespread in high-dimensional statistical inference.

3 Minimax Framework and High-Dimensional Scaling

There are two main ingredients in the minimax framework. The first is the class of probability distributions under consideration. These are usually associated with some parameter space corresponding to the structural constraints. Formally, suppose that $\lambda_{1}>\lambda_{2}$. Then we may write eq. 2 as

The second ingredient in the minimax framework is the loss function. In the case of subspace estimation, an obvious criterion for evaluating the quality of an estimator $\hat{\Theta}$ is the squared distance between $\hat{\Theta}$ and $\Theta$. However, it is not appropriate because $\Theta$ is not unique—$\Theta$ and $\Theta V$ span the same subspace for any $k\times k$ orthogonal matrix $V$. On the other hand, the orthogonal projections $\Pi=\Theta\Theta^{T}$ and $\hat{\Pi}=\hat{\Theta}\hat{\Theta}^{T}$ are unique. So we consider the loss function defined by the Frobenius norm of their difference:

In the case where $k=1$, the only possible non-uniqueness in the leading eigenvector is its sign ambiguity. Still, we prefer to use the above loss function in the form

because it generalizes to the case $k>1$. Moreover, when $k=1$, it turns out to be equivalent to both the Euclidean distance between $\theta_{1}$, $\hat{\theta}_{1}$ (when they belong to the same half-space) and the magnitude of the sine of the angle between $\theta_{1}$, $\hat{\theta}_{1}$. (See Lemmas A.1.1 and A.1.2 in the Appendix.)

Our goal in this work is to provide non-asymptotic bounds on the minimax error

Consider the constrained maximization problem

5 Related Work

Operator norm consistent estimates of the covariance matrix automatically imply consistent estimates of eigenspaces. This follows from matrix perturbation theory [[, see, e.g.,]]StewartAndSun. There has been much work on finding operator norm consistent covariance estimators in high-dimensions under assumptions on the sparsity or bandability of the entries of $\Sigma$ or $\Sigma^{-1}$ [[, see, e.g.,]]Bickel:2008a,Bickel:2008,ElKaroui:2008. Minimax results have been established in that setting by . However, sparsity in the covariance matrix and sparsity in the leading eigenvector are different conditions. There is some overlap (e.g. the spiked covariance model), but in general, one does not imply the other.

In next section, we present our main results along with some additional conditions to guarantee that estimation over $\mathcal{M}_{q}$ remains non-trivial. The main steps of the proofs are in Section 3. In the proofs we state some auxiliary lemmas. They are mainly technical, so we defer their proofs to the Appendix. Section 4 concludes the paper with some comments on extensions of this work.

Main Results

Our minimax results are formulated in terms of non-asymptotic bounds that depend explicitly on $(n,p,R_{q},\lambda_{1},\lambda_{2})$. To facilitate presentation, we introduce the notations

There exists $\alpha\in(0,1]$, depending only on $q$, such that

where $\kappa\leq c\alpha/16$ is a constant depending only on $q$, and

In the high-dimensional case that we are interested, where $p>n$, the condition that

for some $\alpha^{\prime}\in$, is sufficient to ensure that (5) holds for $q\in(0,1]$. Alternatively, if we let $\alpha=1-q/2$ then (5) is satisfied for $q\in(0,1]$ if

The relationship between $n$, $p$, $R_{q}$ and $\sigma^{2}$ described in Assumption 2.1 indicates a regime in which the inference is neither impossible nor trivially easy. We can now state our first main result.

Let $q\in$. If Assumption 2.1 holds, then there exists a universal constant $c>0$ depending only on $q$, such that every estimator $\hat{\theta}_{1}$ satisfies

Random variables with finite $\psi_{\alpha}$-norm correspond to those whose tails are bounded by $\exp(-Cx^{\alpha})$.

The case $\alpha=2$ is important because it corresponds to random variables with sub-Gaussian tails. For example, if $Y\sim\mathcal{N}(0,\sigma^{2})$ then $\lVert Y\rVert_{\psi_{2}}\leq C\sigma$ for some positive constant $C$. See [24, Chapter 2.2] for a complete introduction.

Assumption 2.2 holds for a variety of distributions, including the multivariate Gaussian (with $K^{2}=8/3$) and those of bounded random vectors. Under this assumption, we have the following theorem.

If the distribution of $(X_{1},\ldots,X_{n})$ belongs to $\mathcal{M}_{q}(\lambda_{1},\lambda_{2},\bar{R}_{q},\alpha,\kappa)$ and satisfies Assumptions 2.1 and 2.2, then there exists a constant $c>0$ depending only on $K$ such that the following hold:

Proofs of Main Results

Our main tool for proving the minimax lower bound is the generalized Fano Method . The following version is from [26, Lemma 3].

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

The method works by converting the problem from estimation to testing by discretizing the parameter space, and then applying Fano’s Inequality to the testing problem. (The $\beta_{N}$ term that appears above is an upper bound on the mutual information.)

To be successful, we must find a sufficiently large finite subset of the parameter space such that the points in the subset are $\alpha_{N}$-separated under the loss, yet nearly indistinguishable under the KL divergence of the corresponding probability measures. We will use the subset given by the following lemma.

Fix $\epsilon\in(0,1]$ and let $\Theta_{\epsilon}$ denote the set given by Lemma 3.1.2. With Lemma A.1.2 we have

for all distinct pairs $\theta_{1},\theta_{2}\in\Theta_{\epsilon}$. For each $\theta\in\Theta_{\epsilon}$, let

where $\sigma^{2}=\lambda_{1}\lambda_{2}/(\lambda_{1}-\lambda_{2})^{2}$.

Thus, we have found a subset of the parameter space that conforms to the requirements of Lemma 3.1.1, and so

for all $\epsilon\in(0,1]$. The final step is to choose $\epsilon$ of the correct order. If we can find $\epsilon$ so that

For a constant $C\in(0,1)$ to be chosen later, let

We consider each of the two cases in the above $\min\{\cdots\}$ separately.

Then $\epsilon^{2}=1$ and by rearranging (11)

observe that the function $x\mapsto x\log[(p-1)/x]$ is increasing on $[1,(p-1)/e]$, and, by Assumption 2.1, this interval contains $\bar{R}_{q}^{2/(2-q)}$. If $p$ is large enough so that $p-1\geq\exp\{(4/c)\log 2\}$, then

Thus, eqs. 8 and 9 are satisfied, and we conclude that

as long as $C^{2}\leq c/16$ and $p-1\geq\exp\{(4/c)\log 2\}$.

and it is straightforward to check that Assumption 2.1 implies that if $C^{q}\geq\kappa^{q}$, then there is $\alpha\in(0,1]$, depending only on $q$, such that

where the last inequality is obtained by plugging in (13) and (14).

If we choose $C^{2}\leq c\alpha/16$, then combining (10) and (3.1), we have

and eq. 8 is satisfied. On the other hand, by (12) and the fact that $\bar{R}_{q}\geq 1$, we have

The function $x\mapsto x\log[(p-1)/x^{2/(2-q)}]$ is increasing on $[1,(p-1)^{1-q/2}/e]$ and, by Assumption 2.1, $1\leq\bar{R}_{q}\leq(p-1)^{1-q/2}/e$. If $p-1\geq\exp\{[4/(c\alpha)]\log 2\}$, then

and eq. 9 is satisfied. So we can conclude that

as long as $C^{2}\leq c\alpha/16$ and $p-1\geq\exp\{[4/(c\alpha)]\log 2\}$.

Looking back at cases 1 and 2, we see that because $\alpha\leq 1$, the conditions that $\kappa^{2}\leq C^{2}\leq c\alpha/16$ and $p-1\geq\exp\{[4/(c\alpha)]\log 2\}$ are sufficient to ensure that

for a constant $c^{\prime}>0$ depending only on $q$. ∎

2 Proof of the Upper Bound (Theorem 2.2)

We begin with a lemma that bounds the curvature of the matrix functional $\langle\Sigma,bb^{T}\rangle$.

We consider the cases $q\in(0,1)$, $q=1$, and $q=0$ separately.

By applying Hölder’s Inequality to the right side of eq. 18 and rearranging, we have

Let $t>0$. We can use a standard truncation argument [[, see, e.g.,]Lemma 5]Raskutti:2011 to show that

Letting $t=\lVert\operatorname{vec}(S-\Sigma)\rVert_{\infty}/(\lambda_{1}-\lambda_{2})$ and joining with eq. 19 gives us

If we define $m$ implicitly so that $\epsilon=m\sqrt{2}t^{1-q/2}R_{q}$, then the preceding inequality reduces to $m^{2}/2\leq m+1$. If $m\geq 3$, then this is violated. So we must have $m<3$ and hence

Combining the above discussion with the sub-Gaussian assumption, the next lemma allows us to bound $\lVert\operatorname{vec}(S-\Sigma)\rVert_{\infty}$.

If Assumption 2.2 holds and $\Sigma$ satisfies (2), then there is an absolute constant $c>0$ such that

Combining this with the trivial bound $\epsilon\leq 2$, yields

for an appropriate constant $c>0$, depending only on $K$. This completes the proof for the case $q\in(0,1)$.

2.2 Case 2: q=1𝑞1q=1

The next lemma provides a bound for the supremum.

If Assumption 2.2 holds and $\Sigma$ satisfies (2), then there is an absolute constant $c>0$ such that

Assumption 2.1 guarantees that $R_{1}^{2}\in[1,p/e]$. Thus, we can apply Lemma 3.2.3 and an argument similar to that used with (21) to complete the proof for the case $q=1$.

2.3 Case 3: q=0𝑞0q=0

where $\lVert\,\cdot\,\rVert_{S_{1}}$ denotes the sum of the singular values. Divide both sides by $\epsilon$, rearrange terms, and then take the expectation to get

If Assumption 2.2 holds and $\Sigma$ satisfies (2), then there is an absolute constant $c>0$ such that

Taking $d=2R_{0}$ and applying an argument similar to that used with (21) completes the proof of the $q=0$ case. ∎

Conclusion and Further Extensions

V. Q. Vu was supported by a NSF Mathematical Sciences Postdoctoral Fellowship (DMS-0903120). J. Lei was supported by NSF Grant BCS0941518. We thank the anonymous reviewers for their helpful comments.

References

Appendix A APPENDIX - SUPPLEMENTARY MATERIAL

We state below two results that we use frequently in our proofs. The first is well-known consequence of the CS decomposition. It relates the canonical angles between subspaces to the singular values of products and differences of their corresponding projection matrices.

The singular values of $\Pi_{\mathcal{X}}(I_{p}-\Pi_{\mathcal{Y}})$ are

The singular values of $\Pi_{\mathcal{X}}-\Pi_{\mathcal{Y}}$ are

If in addition $\lVert x-y\rVert_{2}\leq\sqrt{2}$, then

By Lemma A.1.1 and the polarization identity

The upper bound follows immediately. Now if $\lVert x-y\rVert_{2}^{2}\leq 2$, then the above right-hand side is bounded from below by $\lVert x-y\rVert_{2}^{2}/2$. ∎

A.2 Proofs for Theorem 2.1

Our construction is based on a hypercube argument. We require a variation of the Varshamov-Gilbert bound due to . We use a specialization of the version that appears in [15, Lemma 4.10].

Let $d$ be an integer satisfying $1\leq d\leq(p-1)/4$. There exists a subset $\Omega_{d}\subset\{0,1\}^{p-1}$ that satisfies the following properties:

$\lVert\omega\rVert_{0}=d$ for all $\omega\in\Omega_{d}$,

$\lVert\omega-\omega^{\prime}\rVert_{0}>d/2$ for all distinct pairs $\omega,\omega^{\prime}\in\Omega_{d}$, and

$\log\lvert\Omega_{d}\rvert\geq cd\log((p-1)/d)$, where $c\geq 0.233$.

Let $d\in[1,(p-1)/4]$ be an integer, $\Omega_{d}$ be the corresponding subset of $\{0,1\}^{p-1}$ given by preceding lemma,

Clearly, $\Theta$ satisfies the following properties:

$\epsilon/\sqrt{2}<\lVert\theta_{1}-\theta_{2}\rVert_{2}\leq\sqrt{2}\epsilon$ for all distinct pairs $\theta_{1},\theta_{2}\in\Theta_{d}$,

$\lVert\theta\rVert_{q}^{q}\leq 1+\epsilon^{q}d^{(2-q)/2}$ for all $\theta\in\Theta$, and

$\log\lvert\Theta\rvert\geq cd[\log(p-1)-\log d]$, where $c\geq 0.233$.

for all $\theta\in\Theta$. To complete the proof we will show that $\log\lvert\Theta\rvert$ satisfies the lower bound claimed by the lemma. Note that the function $a\mapsto a\log[(p-1)/a]$ is increasing on $[0,(p-1)/e]$ and decreasing on $[(p-1)/e,\infty)$. So if

because $d=\lfloor a\rfloor\geq a/2$. Moreover, since $d\leq(p-1)/4$ and the above right hand side is maximized when $a=(p-1)/e$, the inequality remains valid for all $a\geq 0$ if we replace the constant $(c/2)$ with the constant

Let $A_{i}=x_{i}x_{i}^{T}$ for $i=1,2$. Then $\Sigma_{i}=\lambda_{1}A_{i}+\lambda_{2}(I_{p}-A_{i})$. Since $\Sigma_{1}$ and $\Sigma_{2}$ have the same eigenvalues and hence the same determinant,

The spectral decomposition $\Sigma_{2}=\lambda_{1}A_{2}+\lambda_{2}(I_{p}-A_{2})$ allows us to easily calculate that

Since orthogonal projections are idempotent, i.e. $A_{i}A_{i}=A_{i}$,

Using again the idempotent property and symmetry of projection matrices,

A.3 Proofs for Theorem 2.2

Since $\theta_{1}$ is an eigenvector of $\Sigma$ corresponding to the eigenvalue $\lambda_{1}$,

The last inequality follows from Lemma A.1.1. ∎

Using the elementary inequality $2|ab|\leq a^{2}+b^{2}$, we have by Assumption 2.2 that

In the third line, we used the fact that the $\psi_{1}$-norm is bounded above by a constant times the $\psi_{2}$-norm [[, see]p. 95]vanderVaartAndWellner. By a generalization of Bernstein’s Inequality for the $\psi_{1}$-norm [[, see]Section 2.2]vanderVaartAndWellner , for all $t>0$

This implies [24, Lemma 2.2.10] the bound

Adding LABEL:eq:d-bound and 23 and then adjusting the constant $c$ gives the desired result, because

The result uses ’s generic chaining method, and allows us to reduce the problem to bounding the supremum of a Gaussian process. The statement of the result involves the generic chaining complexity, $\gamma_{2}(B,d)$, of a set $B$ equipped with the metric $d$. We only use a special case, $\gamma_{2}(B,\lVert\,\cdot\,\rVert_{2})$, where the complexity measure is equivalent to the expectation of the supremum of a Gaussian process on $B$. We refer the reader to for a complete introduction.

Let $Z_{i}$, $i=1,\ldots,n$ be i.i.d. random variables. There exists an absolute constant $c$ for which the following holds. If $\mathcal{F}$ is a symmetric class of mean-zero functions then

where $d_{\psi_{1}}=\sup_{f\in\mathcal{F}}\lVert f\rVert_{\psi_{1}}$.

and $D_{2}(b):=\langle\bar{Z},\Sigma^{1/2}b\rangle^{2}$, respectively. To apply Lemma A.3.1 to $D_{1}$, define the class of linear functionals

and we are in the setting of Lemma A.3.1.

First, we bound the $\psi_{1}$-diameter of $\mathcal{F}$.

Next, we bound $\gamma_{2}(\mathcal{F},\psi_{2})$ by showing that the metric induced by the $\psi_{2}$-norm on $\mathcal{F}$ is equivalent to the Euclidean metric on $B$. This will allow us to reduce the problem to bounding the supremum of a Gaussian process. For any $f,g\in\mathcal{F}$, by Assumption 2.2,

where $b_{f},b_{g}\in B$. Thus, by [23, Theorem 1.3.6],

Then applying Talagrand’s Majorizing Measure Theorem [23, Theorem 2.1.1] yields

where we used the assumption that $R_{1}^{2}\leq p/e$ in the last inequality. Now we apply Lemma A.3.1 to get

Turning to $D_{2}(b)$, we can take $n=1$ in Lemma A.3.1 and use a similar argument as above, because

We just need to bound the $\psi_{2}$-norms of $f(\bar{Z})$ and $(f-g)(\bar{Z})$ to get bounds that are analogous to eqs. 24 and 25. Since $\bar{Z}$ is the sum of the independent random variables $Z_{i}/n$,

Putting together the bounds for $D_{1}$ and $D_{2}$ and then adjusting constants completes the proof. ∎

Using a similar argument as in the proof of Lemma 3.2.3 we can show that

and $Y$ is a $p$-dimensional standard Gaussian $Y$. Thus we can reduce the problem to bounding the supremum of a Gaussian process.

Since $\langle Y,b\rangle$ is a standard Gaussian for every $b\in\mathcal{N}$, a union bound [24, Lemma 2.2.2] implies

In the second line, we used the binomial coefficient bound $\binom{p}{d}\leq(ep/d)^{d}$. If we take $\delta=1/4$, then

where we used the assumption that $d<p/2$. Thus,