Sparse PCA: Optimal rates and adaptive estimation

Introduction

Due to dramatic advances in science and technology, high-dimensional data are now routinely collected in a wide range of fields including genomics, signal processing, risk management and portfolio allocation. In many applications, the signal of interest lies in a subspace of much lower dimension and the between-sample variation is determined by a small number of factors. For example, in spectroscopy, the variation of the infrared and ultraviolet spectra is driven by the concentration levels of a small number of chemical components in the system Varmuza09 . In financial econometrics, it is commonly believed that the variation in asset returns is driven by a small number of common factors combined with random noise Chamberlain83 .

Principal component analysis (PCA) is one of the most commonly used techniques in multivariate analysis for dimension reduction and feature extraction, and is particularly well suited for the settings where the data is high-dimensional but the signal has a low-dimensional structure. PCA has a wide array of applications, ranging from image recognition to data compression to clustering. In the conventional setting where the dimension of the data is relatively small compared with the sample size, the principal eigenvectors of the covariance matrix is typically estimated by the leading eigenvectors of the sample covariance matrix which are consistent when the dimension $p$ is fixed, and the sample size $n$ increases Anderson03 . However, in the high-dimensional setting where $p$ can be much larger than $n$, this approach leads to very poor estimates. At various levels of rigor and generality, a series of papers Hoyle04 , Baik06 , Paul07 , Nadler08 , JohnstoneLu09 , Jung09 , Birnbaum12 showed that the sample principal eigenvectors are no longer consistent estimates of their population counterparts. For example, Baik and Silverstein Baik06 and Paul Paul07 showed that if $p/n\to\gamma\in(0,1)$ as $n\to\infty$, and the largest eigenvalue $\lambda_{1}\leq\sqrt{\gamma}$ and is of unit multiplicity, then the leading sample principal eigenvector $\hat{\mathbf{v}}_{1}$ is asymptotically almost surely orthogonal to the leading population eigenvector $\mathbf{v}_{1}$, that is, $|\mathbf{v}_{1}^{\prime}\hat{\mathbf{v}}_{1}|\to 0$ almost surely. Thus, in this case, $\hat{\mathbf{v}}_{1}$ is not useful at all as an estimate of $\mathbf{v}_{1}$. Even when $\lambda_{1}>\sqrt{\gamma}$, the angle between $\mathbf{v}_{1}$ and $\hat{\mathbf{v}}_{1}$ still does not converge to zero unless $\lambda_{1}\to\infty$. In addition to being inconsistent, sample principal eigenvectors have nonzero loadings in all the coordinates. This renders their interpretation difficult when the dimension $p$ is large.

In view of the above negative results in the high-dimensional setting, a natural approach to principal component analysis in high dimensions is to impose certain structural constraint on the leading eigenvectors. One of the most popular assumptions is that the leading eigenvectors have a certain type of sparsity. In this case, the problem is commonly referred to as sparse PCA in the literature. The sparsity constraint reduces the effective number of parameters and facilitates interpretation.

Various regularized estimators of the leading eigenvectors have been proposed in the literature. See, for example, Jolliffe03 , Zou06 , dAspremont07 , ShenHuang08 , Solo08 , Witten09 , Journee10 . Theoretical analysis has so far mainly focused on the rank-one case, that is, estimating the leading principal eigenvector $\mathbf{v}_{1}$. In this case, Johnstone and Lu JohnstoneLu09 showed that the classical PCA performed on a selected subset of variables with the largest sample variances leads to a consistent estimator of $\mathbf{v}_{1}$ if the ordered coefficients of $\mathbf{v}_{1}$ have rapid decay. Shen, Shen and Marron Shen11 and Yuan and Zhang Yuan11 proposed other consistent estimators when $\mathbf{v}_{1}$ has a bounded number of nonzero coefficients. Vu and Lei Vu12 studied the rates of convergence of estimation under various sparsity assumptions on $\mathbf{v}_{1}$, and Lounici Lounici12 further considers the minimax rates with missing data. Amini and Wainwright Amini09 investigated the variable selection property of the methods by JohnstoneLu09 and dAspremont07 when $\mathbf{v}_{1}$ has $k$ nonzero entries all of the same magnitude. Berthet and Rigollet Berthet12 considered minimax detection when $\mathbf{v}_{1}$ has a bounded number of nonzeros.

More recently, for estimating a fixed number $r\geq 1$ of leading eigenvectors as $n,p\to\infty$, Birnbaum et al. Birnbaum12 studied minimax rates of convergence and adaptive estimation of the individual leading eigenvectors when the ordered coefficients of each eigenvector have rapid decay. When $r>1$ and some of the leading eigenvalues have multiplicity great than one, the individual leading eigenvectors can be unidentifiable. On the other hand, the principal subspace spanned by them is always uniquely defined. Ma Ma11 proposed a new method for estimating the principal subspace and derived rates of convergence of the estimator under similar conditions to those in Birnbaum12 .

2 Estimation of principal subspace

In this paper, we focus on the estimation of the principal subspace. Both minimax and adaptive estimation are considered. Throughout the paper, let $\mathbf{X}$ be an $n\times p$ data matrix generated as

Here $\mathbf{U}$ is the $n\times r$ random effects matrix with i.i.d. $N(0,1)$ entries, $\mathbf{D}=\operatorname{diag}(\lambda_{1}^{1/2},\ldots,\lambda_{r}^{1/2})$ with $\lambda_{1}\geq\cdots\geq\lambda_{r}>0$, $\mathbf{V}$ is $p\times r$ orthonormal and $\mathbf{Z}$ has i.i.d. $N(0,\sigma^{2})$ entries which are independent of $\mathbf{U}$. Equivalently, one can think of $\mathbf{X}$ as an $n\times p$ matrix with rows independently drawn from the distribution $N(0,\bolds{\Sigma})$, where the covariance matrix $\bolds{\Sigma}$ is given by

Here $\bolds{\Lambda}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{r})$ and $\mathbf{V}=[\mathbf{v}_{1},\ldots,\mathbf{v}_{r}]$ is $p\times r$ with orthonormal columns. The $r$ largest eigenvalues of $\bolds{\Sigma}$ are $\lambda_{i}+\sigma^{2}$, $i=1,\ldots,r$, and the rest are all equal to $\sigma^{2}$. The $r$ leading eigenvectors of $\bolds{\Sigma}$ are given by the columns of $\mathbf{V}$. Since the spectrum of $\bolds{\Sigma}$ has $r$ spikes, the covariance structure (2) is commonly known as the spiked covariance matrix model Johnstone01 in the literature.

which is a commonly used metric to gauge the distance between linear subspaces. It also coincides with twice the sum of the squared sines of the principal angles between the respective linear span.

In particular, if $q=0$, then we have $1\leq r\leq s\leq p$. Throughout the paper, we assume that (8) holds.

3 Optimal rates of convergence

Combining the upper and lower bound results developed in Section 2, we establish the following minimax rates of convergence for estimating the principal subspace $\operatorname{span}(\mathbf{V})$ under the loss (3). We focus here on the exact sparse case of $q=0$; the optimal rates for the general case of $q\in(0,2)$ are given in Section 2. For two sequences of positive numbers $a_{n}$ and $b_{n}$, we write $a_{n}\gtrsim b_{n}$ when $a_{n}\geq cb_{n}$ for some absolute constant $c>0$ and $a_{n}\lesssim b_{n}$ when $b_{n}\gtrsim a_{n}$. Finally, we write $a_{n}\asymp b_{n}$ when both $a_{n}\gtrsim b_{n}$ and $a_{n}\lesssim b_{n}$ hold.

as long as the right-hand side of (9) does not exceed some absolute constant. Otherwise, there exists no consistent estimator.

The rate of convergence in (9) depends optimally on all the parameters $s,p,r,n$ and $\lambda$. The result thus provides a precise characterization of the difficulty of the principal subspace estimation problem in terms of the minimax rates over a wide range of parameter values.

We then construct an explicit estimator using an aggregation scheme, which is shown to attain the same rates of convergence as those of the minimax lower bounds. The matching lower and upper bounds together establish the optimal rates of convergence. This aggregation method can potentially be useful for other high-dimensional sparse PCA problems as well. Aggregation methods have been widely used and well studied in statistics literature. See, for example, Juditsky and Nemirovski JN00 , Yang Yang00 , Nemirovski Nemirovski00 and Rigollet and Tsybakov RT11 . To the best of our knowledge, this is the first application of the aggregation approach to sparse PCA which yields optimality results.

4 Adaptive estimation

The rate-optimal aggregation estimator depends on the model parameters that are usually unknown in practice and is unfortunately not computationally feasible when $p$ is large. We then propose an adaptive estimation procedure that is fully data driven and easily implementable. The estimator is shown to attain the optimal rate of convergence simultaneously over a large collection of the parameter spaces defined in (1.2).

The proposed method is based on a reduction scheme. By a conditioning argument, the original sparse PCA problem is reduced to a high-dimensional regression problem with orthogonal design and group sparsity on the regression coefficients. Then, we apply the model selection penalty idea from Birge01 to construct the final estimator.

A key step in the reduction scheme is the construction of two new samples in the form of (1), which share the same realization of the random effects $\mathbf{U}$ but have independent copies of the noise matrix $\mathbf{Z}$. This construction works because a common realization of $\mathbf{U}$ is critical in guaranteeing a sufficient signal-to-noise ratio in the resulting regression problem. In contrast, splitting the original sample into two halves fails to achieve this goal. On the other hand, the independence of the noise components ensures that the regression problem has white noise structure. The adaptivity and minimax optimality of the subspace estimator depend heavily on those of the regression coefficient estimator. Thus, as a byproduct of the analysis, we also show that our estimator for regression coefficients is adaptively rate optimal under group sparsity. To the best of our knowledge, the specific estimator and its adaptive optimality is also new in the literature.

5 Other related work

The present paper is related to a fast growing literature on estimating sparse covariance/precision matrices as well as low-rank matrices. Significant advances have been made on optimal estimation of the whole covariance or precision matrix. Many regularization methods, including banding, tapering, thresholding and penalization, have been proposed. In particular, Cai, Zhang and Zhou CZZ10 established the optimal rate of convergence for estimating a class of bandable covariance matrices under the spectral norm. Cai and Yuan CY12 proposed a block thresholding procedure which is shown to be adaptively rate-optimal over a wide range of collections of bandable covariance matrices. Bickel and Levina BJ08b introduced a thresholding procedure and obtained rates of convergence for sparse covariance matrix estimation. Cai and Zhou CZ12 established the minimax rates of convergence for estimating sparse covariance matrices under a range of matrix norms including the spectral norm. Cai, Liu and Zhou CLZ12 obtained the optimal rate of convergence for estimating the sparse precision matrices.

Our work is also related to another active area of research, namely, the recovery of low-rank matrices based on noisy observations. Negahban and Wainwright Negahban11 studied (near) low-rank matrix recovery by $M$-estimators under restricted strong convexity based on the penalized nuclear norm minimization over matrices. Koltchinskii, Lounici and Tsybakov Koltchinskii11 considered estimation of low-rank matrices based on a trace regression model which includes matrix completion as a special case. A nuclear norm penalized estimator was proposed and a general sharp oracle inequality was established. See also Recht, Fazel and Parrilo rankmin and Rohde and Tsybakov Rhode11 .

6 Organization of the paper

The rest of the paper is organized as follows. After introducing basic notation, Section 2 establishes the minimax rates of convergence for estimating the principal subspace by obtaining matching minimax lower and upper bounds. An aggregation estimator is constructed and shown to be rate optimal. Section 3 introduces an adaptive estimation procedure for the principal subspace which is fully data driven and easily computable. It is shown that this estimator attains the optimal rates of convergence simultaneously over a large collection of parameter spaces. Connections to other related problems are discussed in Section 5. The proofs of the main results and key technical lemmas are given in Section 6 and some additional technical arguments are contained in the supplementary material appsm .

Minimax rates for principal subspace estimation

We establish in this section the minimax rates of convergence for estimating the principal subspace in two steps. First, minimax lower bounds are obtained for the estimation problem under the loss (3). Then an aggregation estimator is introduced and is shown to attain the same rates as given in the lower bounds, under mild conditions on the parameters. The matching lower and upper bounds thus establish the minimax rates of convergence.

We first establish the minimax lower bounds which are instrumental in obtaining the optimal rates of convergence. In view of the upper bounds to be given in Section 2.2 by an aggregation procedure, these lower bounds are minimax rate optimal under mild conditions.

Before proceeding to the precise statements, we introduce the following notation: let

$k_{q}^{*}=p$ if and only if $s\geq p(\frac{r+1}{nh(\lambda)})^{q/2}$, in which case the effective dimension coincides with the ambient dimension.

$s\mapsto k_{q}^{*}$ is increasing. Moreover, there exists a function $\tau_{q}$, such that $k^{*}_{q}(as,p,r,n,\lambda)\leq k^{*}_{q}(s,p,r,n,\lambda)\tau_{q}(a)$ for any $a\geq 1$.

$k_{q}^{*}\gtrsim s$ if and only if the assumption (16) holds. See Figure 1 for a graphical illustration on the dependence of the effective dimension $k_{q}^{*}$ on various parameters.

Without loss of generality, we assume unit noise variance ($\sigma^{2}=1$) from now on. All results hold for a general $\sigma$ by replacing $\lambda$ with $\lambda/\sigma^{2}$. We consider the lower bounds separately in two cases: $0<q<2$ and $q=0$.

for some sufficiently large absolute constant $C_{0}$. Then there exists a constant $c$ depending only on $q$ and an absolute constant $c_{0}$, such that the minimax risk for estimating $\mathbf{V}$ over the parameter space $\Theta=\Theta_{q}(s,p,r,\lambda)$ satisfies

For the case of $q=0$ we have the following lower bound:

Let $p,s,r$ be integers such that $1\leq r\leq s\leq p$. Let the observed matrix $\mathbf{X}$ be generated by model (1) with $\sigma=1$. Then the minimax risk for estimating $\mathbf{V}$ over the parameter space $\Theta=\Theta_{0}(s,p,r,\lambda)$ satisfies

2 Optimal estimation via aggregation

We now show that the lower bounds given in Section 2.1 are indeed rate optimal under mild technical conditions. The optimal estimator of $\mathbf{V}$ is constructed using sample splitting and aggregation. The estimator is theoretically interesting but computationally intensive. We will construct a data-driven and easily implementable estimator in Section 3 under stronger conditions.

We first note that the loss function (3) satisfies

Moreover, the loss function is invariant under orthogonal complement, that is, $L(\mathbf{V},\widehat{\mathbf{V}})=L(\mathbf{V}^{\perp},\widehat{\mathbf{V}}^{\perp})$, where $[\mathbf{V},\mathbf{V}^{\perp}],[\widehat{\mathbf{V}},\widehat{\mathbf{V}}^{\perp}]$ are orthogonal matrices. Therefore the loss (19) admits the following upper bound:

For notational simplicity we assume that the sample size is $2n$ and we split the sample equally according to \mathbf{X}=\bigl{[}{{\matrix{\mathbf{X}_{(1)}\cr\mathbf{X}_{(2)}}}}\bigr{]}, where $\mathbf{X}_{(i)}=\mathbf{U}_{(i)}\mathbf{D}\mathbf{V}^{\prime}+\mathbf{Z}_{(i)},i=1,2$. Denote by $\mathbf{S}_{(i)}=\frac{1}{n}\mathbf{X}_{(i)}^{\prime}\mathbf{X}_{(i)}$ the corresponding sample covariance matrix. The main idea is to construct a family of estimators $\{\widehat{\mathbf{V}}_{B}\}$ using the first sample, indexed by the row support $B\subset[p]$, where $\widehat{\mathbf{V}}_{B}$ is the optimal estimator one would use if one knew beforehand that $\operatorname{supp}(\mathbf{V})=B$. Then we aggregate these estimators by selection using the second sample.

Recall the effective dimension $k^{*}_{q}$ defined in (13). For each $B\subset[p]$ such that $|B|=k^{*}_{q}$, we define $\widehat{\mathbf{V}}_{B}\in O(p,r)$ as the $r$ leading singular vectors of $\mathbf{J}_{B}\mathbf{S}_{(1)}\mathbf{J}_{B}$, where $\mathbf{J}_{B}$ is the diagonal matrix given by

Given the collection of the $\widehat{\mathbf{V}}_{B}$’s, we set

It is natural to use the same sample covariance matrix to construct the $\widehat{\mathbf{V}}_{B}$’s and to select $B^{*}$. The main advantage of sample splitting is to decouple the selection of the support and the computation of the estimator. Thus, conditioning on the first sample, we can treat the candidate estimators as if they are deterministic, which greatly facilitates the analysis. Sample splitting is commonly used in aggregation based estimation, where a sequence of estimators is constructed from the first sample and the second sample is used to aggregate these candidates to produce a final estimator.

Let $q\in[0,2)$. Let $k^{*}_{q}$ be defined in (13). Let $\widehat{\mathbf{V}}_{*}$ be the aggregated estimator defined in (23). Assume that

for some sufficiently large constant $C_{0}$. Then there exists a constant $C$ depending only on $\kappa$ and $q$ such that for $\Theta=\Theta_{q}(s,r,p,\lambda)$,

where $\Psi(k,p,r,n,\lambda)$ and $k^{*}_{q}$ are defined in (11) and (13), respectively. Moreover, if $q=0$, then $\Psi$ in (27) can be replaced by $\Psi_{0}$ defined in (12) with $k^{*}_{0}=s$, and condition (25) can be dropped.

When $q\in(0,2)$, under the conditions of Theorems 2 and 4, the lower and upper bounds together yield the minimax rates of convergence $\Psi(k^{*}_{q},p,r,n,\lambda)$ given in (11) with the optimal dependence on all the parameters, in particular the eigenvalues and the rank. When $q=0$, the lower and upper bounds match under less restrictive conditions, which will be discussed in more detail in Remark 2 below.

3 Comments

We conclude this section with a few important remarks.

Comparing the lower and upper bounds for $q=0$ in Theorems 2 and 4, we see a sufficient condition for the minimax rate to match (and hence coincide with $\Psi_{0}$) is

It is interesting to note that under the condition (28), the minimax rate for estimating the $r$ leading singular vectors depend on the $r$ only through $r(s-r)$, which is the dimension of the Grassmannian manifold $G(s,r)$. Therefore the dependence on $r$ is not monotonic, with the worst case happening at $r=\frac{s}{2}$. However, it should be noted that in order for the minimax rate to coincide with $\Psi_{0}$, it is necessary to have $r$ strictly bounded away from $s$, for example, in the regime of (28). When $r=s$, the lower bound in Theorem 3 becomes zero. In this degenerate case, the only uncertainty is in the support of $\mathbf{V}$. The minimax rate is indeed much faster than $\Psi_{0}$, because in this regime the support can be estimated accurately. See Section 7.3 in the supplementary material appsm .

For $q\in(0,2)$, the minimax rate $\Psi$ depends on the effective dimension $k_{q}^{*}$ which is defined implicitly through equations (13)–(14). It is possible to obtain an explicit formula of the minimax rate in some regime. For example, if $s\geq p^{1-\epsilon}(\frac{r+\log p}{nh(\lambda)})^{q/2}$ for some constant $\epsilon\in(0,1)$, then the effective dimension satisfies $k_{q}^{*}\leq p^{1-\epsilon}$. Moreover, we have $k_{q}^{*}\asymp s(\frac{nh(\lambda)}{r+\log p})^{q/2}$. Hence the minimax rate is given by

An interesting side product of the proofs of Theorems 3 and 4 is the following nonasymptotic minimax rate for the regular PCA problem without structural assumptions on the principle subspaces. It is a classical result (see, e.g., Stein56 , Eaton70 ) that when $p\leq n$, the sample covariance matrix is not exact minimax optimal for estimating the whole covariance matrix under certain losses (e.g., the Stein loss). As shown in the next theorem, in the unstructured case, it turns out that the sample version of the principle subspace is minimax rate optimal even in high dimensions. For more details see Theorems 8 and 9 in Sections 6.1 and 6.2.

Let $\Theta=\Theta_{0}(p,p,r,\lambda)$. Let $n\geq C_{0}(r+\log\lambda)$ and $\lambda\geq\break C_{0}\sqrt{{\log(n)}/{n}}$ for some sufficiently large constant $C_{0}$. Then for all $r\in[p]$,

which can be attained by $\widehat{\mathbf{V}}$ consisting of the $r$ leading eigenvectors of the sample covariance matrix $\mathbf{S}$.

Theorem 5 implies that, without structural assumptions on the principle subspace $\mathbf{V}$, consistent estimators exist if and only $\frac{nh(\lambda)}{r(p-r)}\to\infty$. Moreover, unless $nh(\lambda)$ exceeds a constant factor of $p$, even the optimal estimator is within a constant factor of $r\wedge(p-r)$, the upper bound of the loss function.

In the special case of $r=1$, a similar combinatorial procedure to (22)–(23) has been proposed in Vu12 . Using Mendelson’s results on empirical processes Mendelson10 , this procedure requires no sample splitting but can only be shown to attain a convergence rate that is suboptimal in $\lambda$ Vu12 , Theorem 2.2: with $\lambda\to\infty$ and all the other parameters fixed, the upper bound in Vu12 does not vanish. In contrast, the optimal rate $\Psi$ decays at the rate $\lambda^{-(1-q/2)}$ when $k_{q}^{*}<p$ and $\lambda^{-1}$ when $k_{q}^{*}=p$. Comparing with the analysis in Vu12 , the proof of Theorem 4 is more elementary. By exploring the structure of the difference between the sample covariance matrix and the true covariance matrix, we obtained an upper bound that is optimal in all parameters.

Adaptive estimation

The aggregation estimator constructed in Section 2.2 has been shown to be rate optimal. However, it depends on the unknown parameters and is computationally infeasible when $p$ is large. We construct in this section an adaptive estimation procedure for principal subspaces which is fully data driven and easily computable. Furthermore, it is shown that the estimator attains the optimal rate of convergence simultaneously over a large collection of the parameter spaces defined in (1.2).

Let $\mathbf{S}^{i}=\frac{1}{n}(\mathbf{X}^{i})^{\prime}\mathbf{X}^{i}$, $i=0,1$, be the sample covariance matrices for the two samples.

We use the sample $\mathbf{X}^{0}$ to compute an initial estimator $\mathbf{V}^{0}$. A specific procedure for computing the initial estimator $\mathbf{V}^{0}$ will be given in Section 3.2.

where $\mathbf{Y}=\frac{1}{\sqrt{2}}(\mathbf{X}^{1})^{\prime}\mathbf{X}^{0}\mathbf{V}^{0}\mathbf{R}\mathbf{C}^{-1}$, $\bolds{\Theta}=\frac{1}{\sqrt{2}}\mathbf{V}\mathbf{A}\mathbf{R}\mathbf{C}^{-1}$ and $\mathbf{E}=\frac{1}{\sqrt{2}}(\mathbf{Z}^{1})^{\prime}\mathbf{L}$. We shall treat (32) as a regression problem, where $\mathbf{Y}$ is the observed matrix, $\bolds{\Theta}$ is the signal matrix of interest and $\mathbf{E}$ is the additive noise matrix. Equivalently, we think of $\bolds{\Theta}$ as the coefficient matrix, and the design matrix is $\mathbf{X}=\mathbf{I}_{p}$. The reason why this is plausible will be detailed in Section 3.2.

Given $\mathbf{Y}$, we propose the following method for computing $\widehat{\bolds{\Theta}}$. Define

Fix an arbitrary $\delta\in(0,1)$. With slight abuse of notation, define

Then the estimator for $\bolds{\Theta}$ is defined as

Such a penalized least squares approach has been widely used in orthogonal regression with various choices of the penalty functions. See, for example, Birgé and Massart Birge01 and Abramovich et al. ABDJ06 .

In case of multiple minimizers, $\hat{k}$ is chosen to be the smallest one. It is also clear that $\hat{k}$ is easy to compute. With $\hat{k}$, the estimator $\widehat{\bolds{\Theta}}$ is given by $\widehat{\bolds{\Theta}}=[\hat{\bolds{\theta}}_{1},\ldots,\hat{\bolds{\theta}}_{p}]^{\prime}$ where

Note that $\hat{k}$ can be equivalently defined as $\operatorname{argmin}_{k\in[p]}\sum_{i=1}^{k}[(1+\delta)^{2}t_{i}-\|\mathbf{y}_{(i)}\|^{2}]$. Therefore $\|\mathbf{y}_{(\hat{k})}\|^{2}>(1+\delta)^{2}t_{\hat{k}}$ and $\|\mathbf{y}_{(\hat{k}+1)}\|^{2}\leq(1+\delta)^{2}t_{\hat{k}+1}$. Since $t_{k}$ is strictly decreasing in $k$, we obtain that $\|\mathbf{y}_{(1)}\|^{2}\geq\cdots\geq\|\mathbf{y}_{(\hat{k})}\|^{2}>(1+\delta)^{2}t_{\hat{k}}\geq\|\mathbf{y}_{(\hat{k}+1)}\|^{2}\geq\cdots.$ Thus, $|\operatorname{supp}(\widehat{\bolds{\Theta}})|=\hat{k}$.

Step 4: Final estimation. Last but not least, we obtain the estimator $\widehat{\mathbf{V}}$ for $\mathbf{V}$ by orthonormalizing the columns of $\widehat{\bolds{\Theta}}$. The orthonormalization can be completed by the Gram–Schmidt procedure or QR factorization. The estimated subspace is $\operatorname{span}(\widehat{\mathbf{V}})=\operatorname{span}(\widehat{\bolds{\Theta}})$.

An important feature of the above reduction scheme is that the two samples $\mathbf{X}^{0}$ and $\mathbf{X}^{1}$ share the same realization of random factors $\mathbf{U}$ and their only difference is in the noise matrices $\mathbf{Z}^{0}$ and $\mathbf{Z}^{1}$. This is critical for maintaining the right level of signal-to-noise ratio in the regression problem (32). In contrast, splitting the original sample into two halves as in Section 2.2 does not achieve this goal here. Since our analysis relies on the independence of $\mathbf{Z}^{0}$ and $\mathbf{Z}^{1}$, the normality of the noise is crucial to this scheme.

2 Sparse PCA and regression with group sparsity

We now apply the general reduction scheme to the principal subspace estimation problem with the parameter spaces defined in (1.2). In what follows, we first introduce and study a specific estimator for the initial estimation step. Then we derive properties of the proposed estimator for the regression with group sparsity problem. Furthermore, we show that the general reduction scheme paired with the two specific estimators leads to a final estimator which adaptively achieves the optimal rates of estimation over a large collection of the parameter spaces of interest. For clarity of exposition, we regard the rank $r$ as given when introducing the estimators. Data-driven choice of $r$ is discussed at the end of this subsection.

Let $p_{n}\triangleq p\vee n$. We construct the initial estimator $\mathbf{V}^{0}$ via the diagonal thresholding method JohnstoneLu09 as follows:

where $\{s^{0}_{jj}\}_{j=1}^{p}$ are the diagonal elements of $\mathbf{S}^{0}=\frac{1}{n}(\mathbf{X}^{0})^{\prime}\mathbf{X}^{0}$, and $\alpha>0$ is a tuning parameter.

Compute the first $r$ eigenvectors $\{\hat{\mathbf{v}}^{J}_{1},\ldots,\hat{\mathbf{v}}^{J}_{r}\}$ of the submatrix $\mathbf{S}^{0}_{JJ}$.

The following result, proved in Section 7.5 in the supplementary material appsm , gives sufficient conditions on the model parameters and the choice of $\alpha$ to guarantee that the initial estimator $\mathbf{V}^{0}$ is reasonably close to $\mathbf{V}$, which suffices for the initialization of our scheme.

Suppose that $\log n\geq M_{0}\log\lambda$ for some constant $M_{0}>0$. Suppose that

for a sufficiently large constant $C_{0}>0$. If $\mathbf{V}^{0}$ is defined in (38) with a sufficiently large $\alpha\geq\sqrt{10(1+1/M_{0})}$ in (37), then uniformly over $\Theta=\Theta_{q}(s,p,\break r,\lambda)$, we have

hold with probability at least $1-C/[nh(\lambda)]$, where $k_{q}^{*}$ is defined in (13).

We note that condition (40) is critical in establishing the second claim in (41), which ensures that $\mathbf{V}^{0}$ is a reasonable estimator of $\mathbf{V}$. Such a condition is needed for diagonal thresholding to work even when $r=1$. See, for example, condition C3 in Paul05 , page 95. Theorem 4.1 of Birnbaum12 showed that diagonal thresholding could be suboptimal even under a stronger condition than (40).

When $M_{0}$ in Proposition 1 is unknown, we replace it by

where $\sigma_{1}(\mathbf{S}^{0})$ is the largest eigenvalue of $\mathbf{S}^{0}$. This estimate works because $\sigma_{1}(\mathbf{S}^{0})-2$ is an over-estimate of $\lambda$ with high probability Paul07 , Nadler08 , since the noise variance here is two. The estimator (42) allows us to choose $\alpha$ in (37) without explicit knowledge of $M_{0}$.

Orthogonal regression with group sparsity

We first explain why we can treat (32) as a regression problem. When we condition on the values of $\mathbf{U}$ and $\mathbf{Z}^{0}$, the matrix $\mathbf{X}^{0}$ becomes deterministic. Thus, as deterministic functions of $\mathbf{X}^{0}$, the matrices $\mathbf{V}^{0},\mathbf{B},\mathbf{L},\mathbf{C}$ and $\mathbf{R}$ are also deterministic. Furthermore, $\mathbf{A}$ and hence $\bolds{\Theta}$, as deterministic functions of $\mathbf{U}$ and $\mathbf{B}$, are also deterministic. On the other hand, $\mathbf{Z}^{1}$ is independent of both $\mathbf{U}$ and $\mathbf{Z}^{0}$ and hence is independent of $\mathbf{X}^{0}$, $\mathbf{B}$ and $\mathbf{L}$. Thus, the conditional distribution of $\mathbf{Z}^{1}$ on $(\mathbf{U},\mathbf{Z}^{0})$ always has i.i.d. $N(0,2)$ entries, and so the conditional distribution of $\mathbf{E}$ has i.i.d. standard normal entries. Therefore, when we condition on the values of $\mathbf{U}$ and $\mathbf{Z}^{0}$, problem (32) indeed reduces to a standard multivariate regression problem with orthogonal design and white noise.

When the sparsity of $\mathbf{V}$ is specified as in (1.2), we need to consider the following parameter space for $\bolds{\Theta}$:

with $q\in[0,2)$. The parameter $s^{\prime}$ is typically different from $s$ in (1.2), as it also depends on the other model parameters as well as the realization of $\mathbf{U}$ and $\mathbf{Z}^{0}$. However, this will not cause any difficulty in practice, because the estimator proposed in (35) and the associated theorem below remain valid for all values of $s^{\prime}>0$. In the literature of high-dimensional regression, (43) is usually referred to as the group sparsity constraint on the regression coefficients $\bolds{\Theta}$.

For the estimator $\widehat{\bolds{\Theta}}$ in (35), we have following upper bound on its risk. By the lower bounds in Lounici11 for $q=0$, the rates in Theorem 6 are optimal.

for $t_{k}$ defined in (33), and if the set in (44) is empty, we set $k^{\prime}=p$.

Adaptation

With the above preparation, we are now ready to show that if we start with a proper initial estimator $\mathbf{V}^{0}$ [such as that in (38)] and estimate $\bolds{\Theta}$ by (35), then the estimator $\widehat{\mathbf{V}}$ resulting from orthonormalizing the columns of $\widehat{\bolds{\Theta}}$ achieves the optimal rates of convergence. We state the theorem in a slightly more general format. In particular, it holds for the initial estimator in (38) under the conditions of Proposition 1.

Let $\lambda\geq C_{0}$ for some sufficiently large constant $C_{0}$. Let $\Theta=\Theta_{q}(s,p,r,\lambda)$ satisfy the conditions in Theorem 4. Suppose that there exists an initial estimator $\mathbf{V}^{0}$ which satisfies (41) with probability at least $1-C^{\prime}/(nh(\lambda))$. Then the estimator $\widehat{\mathbf{V}}$ obtained by orthonormalizing $\widehat{\bolds{\Theta}}$ in (35) with $\beta>2$ in (33) and $\delta\in(0,1)$ in (34) satisfies

where $k_{q}^{*}$ is defined in (13), and $C>0$ is a constant depending only on $q,\beta$ and $\delta$.

We note that the assumption $\lambda>C_{0}$ is imposed to ensure that the “whitening” procedure in step 3 of the reduction scheme can be performed.

It is interesting to compare the statement of Theorem 7 to the minimax lower bound in Theorems 2–3 as well as the performance of the combinatorial aggregation estimator $\widehat{\mathbf{V}}^{*}$ established in Theorem 4. For any parameter space $\Theta=\Theta_{q}(s,p,r,\lambda)$ such that the conditions of Proposition 1 hold, we could use the $\mathbf{V}^{0}$ in (38), and the resulting $\widehat{\mathbf{V}}$ is guaranteed to achieve the optimal rates of convergence on $\Theta$, which matches the performance of the aggregation estimator for any $q>0$. Moreover, in this case both $\mathbf{V}^{0}$ and $\widehat{\mathbf{V}}$ can be efficiently computed. Hence $\widehat{\mathbf{V}}$ can be used in practice while $\widehat{\mathbf{V}}^{*}$ is computationally intensive. However, in the exact sparse case of $q=0$, the upper bound in Theorem 7 depends on the rank $r$ linearly through $sr$, while the true minimax rate in Theorem 3 depends on $r$ quadratically through $r(s-r)$, which is smaller than $rs$ if $s-r$ is small. The suboptimality of $\widehat{\mathbf{V}}$ in this specific regime is partially due to the fact that our reduction scheme transforms the problem into a regression problem without taking account of the orthogonality structure of the parameter space.

Theorem 7 shows that any estimator $\mathbf{V}^{0}$ satisfying (41) can be used to produce an adaptive estimator. So the task of constructing adaptive optimal estimators is reduced to constructing a “reasonable” estimator.

Consistent estimator of r𝑟r

Last but not least, we discuss how to construct a consistent estimator of $r$ based on data. To this end, recall the definition of the set $J$ in (37), and the matrix $\mathbf{S}^{0}_{JJ}$. We propose to estimate $r$ by

where for any $m>0$ and $M_{0}$ in the conditions of Proposition 1, we define

For this estimator, we have the following result.

Under the condition of Proposition 1, $\hat{r}=r$ holds with probability at least $1-C[nh(\lambda)]^{-1}$.

Under the conditions of Proposition 1 and Theorem 7, Proposition 2 implies that the conclusion in Theorem 7 still holds if we replace $r$ by $\hat{r}$.

Numerical experiments

In this section, we report simulation results comparing the adaptive method proposed in Section 3 with the iterative thresholding method proposed in Ma Ma11 .

In all the results reported here, the sample size $n=1000$ and the ambient dimension $p=2000$. We focus on the case of exact sparsity, that is, $q=0$. The sparsity parameter $s$ takes value in $\{40,80,120,160,200\}$, and the rank $r$ takes value in $\{1,5,10,20\}$. For each $(s,r)$ combination, the $\mathbf{V}$ matrix is obtained from orthonormalizing an $p\times r$ matrix $\mathbf{M}$ where $\mathbf{M}_{i*}$ have i.i.d. $N(0,i^{4})$ entries for $i=1,\ldots,s$ and $\mathbf{M}_{i*}=0$ for all $i>s$. We set the variances of different rows to be different so that the ordered norms of the nonzero rows in $\mathbf{V}$ also exhibit fast decay. When $r=1$, the spike size $\lambda_{1}=20$. When $r>1$, the $\lambda_{i}$’s take $r$ equispaced values such that $\lambda_{r}=10$ and $\lambda_{1}=20$.

When implementing the method in Section 3, we take $\alpha=3$ in (37), $\beta=2.1$ in (33) and $\delta=0.05$ in (34) in all the simulations reported here. In addition, we made a slight modification to the proposed method to obtain better numerical results. We first run the method to obtain an estimator, denoted by $\widehat{\mathbf{V}}_{1}$. Then we switch the roles of $\mathbf{X}^{0}$ and $\mathbf{X}^{1}$ and run the proposed procedure again to obtain a second estimator $\widehat{\mathbf{V}}_{2}$. Finally, we use the $r$ leading eigenvectors of $\widehat{\mathbf{V}}_{1}\widehat{\mathbf{V}}_{1}^{\prime}+\widehat{\mathbf{V}}_{2}\widehat{\mathbf{V}}_{2}^{\prime}$ as the columns of the final estimator $\widehat{\mathbf{V}}$. By Theorem 10 in Section 7.11 in the supplementary material appsm , we have

Here, the first inequality holds because $\sigma_{r}(\widehat{\mathbf{V}}_{1}\widehat{\mathbf{V}}_{1}^{\prime}+\widehat{\mathbf{V}}_{2}\widehat{\mathbf{V}}_{2}^{\prime})\geq\sigma_{r}(\widehat{\mathbf{V}}_{1}\widehat{\mathbf{V}}_{1}^{\prime})=1$ and $\sigma_{r+1}(2\mathbf{V}\mathbf{V}^{\prime})=0$, while the second is by the triangle inequality. By the last display, the theoretical results in Section 3, which apply to both $\widehat{\mathbf{V}}_{1}$ and $\widehat{\mathbf{V}}_{2}$, also apply to the final estimator $\widehat{\mathbf{V}}$. When implementing the iterative thresholding method in Ma Ma11 , we set all tuning parameters at their recommended values.

Table 1 summarizes the average squared Frobenius losses of the proposed method (RegSPCA) and the iterative thresholding method (ITSPCA) over $50$ repetitions for each $(s,r)$ combination. Table 1 shows that for all values of the sparsity parameter, RegSPCA outperformed ITSPCA when $r=5,10$ or $20$, while ITSPCA led to smaller average losses when $r=1$. This demonstrates the competitiveness of RegSPCA in the group sparse setting considered in the present paper. On the other hand, we note that ITSPCA was not designed specifically for handling the group sparsity structure which is the case when $r>1$, and hence its underperformance is not unexpected.

Discussions

We have focused in the present paper on the estimation of the principal subspace $\operatorname{span}(\mathbf{V})$ under the loss (3). The minimax rates of convergence are established and a computationally efficient adaptive estimator is constructed.

Both the current paper and Ma Ma11 consider the problem of sparse subspace estimation under the spiked model, but they differ in several important ways. First, in addition to the sparsity constraint on the leading eigenvectors, the current paper requires them to share support. This extra assumption is motivated by real data applications. For instance, if the observed vectors are the leading Fourier coefficients of random functions with a common covariance kernel, then we expect the leading eigenvectors to have large coefficients only at low frequency coordinates so that the resulting leading eigen-functions in the time domain are smooth. Second, Ma Ma11 focused on the error upper bounds of an adaptive estimator with the subspace rank $r$ assumed to be a fixed constant. Whether the dependence of the bounds on $r$ is optimal was not studied. The current paper conducts an investigation on the dependence of the minimax rates on key model parameters, including $r$ which can grow with $n$ and $p$. Last but not least, we have focused exclusively on the subspace $\operatorname{span}(\mathbf{V})$ which is natural when the spikes are of the same order, while Ma Ma11 considered estimating subspaces spanned by the first few rather than all columns of $\mathbf{V}$. The optimal rates of the latter estimation problem is of most interest when the spikes scale at different rates with $n$ and $p$, which we leave as an interesting problem for future research.

A problem closely related to principal subspace estimation is the estimation of the whole covariance matrix $\bolds{\Sigma}$ under the same structural assumption (1.2). Both minimax estimation and adaptive estimation are of significant interest. Results on minimax rates under the spectral norm loss $L(\widehat{\bolds{\Sigma}},\bolds{\Sigma})=\|\widehat{\bolds{\Sigma}}-\bolds{\Sigma}\|^{2}$ can be found in CMW13 .

It should be noted that our analysis in this paper relies on the normality of the model, which allows us to express the sample in the form of (1). In particular the adaptive procedure requires the independence of $\mathbf{Z}^{0}$ and $\mathbf{Z}^{1}$, which is a consequence of the normality of the noise. It is unclear whether the same results hold for all noise distributions with sub-Gaussian tails. It is an interesting problem to study the robustness of the adaptive procedure and to extend the results to other noise distributions.

Proofs

In this section we prove Theorems 3, 4 and 7. The proofs of the other results, together with those of the key lemmas and some additional technical arguments, are given in the supplementary material appsm .

We first give a lower bound on the oracle risk where we know beforehand the row support of $\mathbf{V}$. This corresponds to a $k$-dimensional unstructured PCA problem, where the goal is to estimate the $r$ leading singular vectors of the covariance matrix. In view of the upper bound in Theorem 9, the rates are minimax optimal.

Let $\Theta=\Theta_{0}(k,k,r,\lambda)$. Then

To prove Theorem 8, we use a minimax lower bound due to Yang and Barron YB99 , Section 7, via local metric entropy, which in turn relies on an argument by Birgé Birge83 . For completeness, we state the result in Proposition 3 and provide a short proof in Section 7.8 in the supplementary material appsm . The method of local metric entropy in an $\frac{1}{\sqrt{n}}$-neighborhood dates back to Le Cam LeCam73 . The advantage of this method is that it only relies on the analytical behavior of the metric entropy of the parameter space, thus allowing us to sidestep constructing explicit packing set in the parameter space.

Let $(\Theta,\rho)$ be a totally bounded metric space and $\{P_{\theta}\dvtx\break\theta\in\Theta\}$ a collection of probability measures. For any $E\subset\Theta$, denote by $\mathcal{N}(E,\epsilon)$ the $\epsilon$-covering number of $E$, that is, the minimal number of balls of radius $\epsilon$ whose union contains $E$. Denote by $\mathcal{M}(E,\epsilon)$ the $\epsilon$-packing number of $E$, that is, the maximal number of points in $E$ whose pairwise distance is at least $\epsilon$. Put

If there exist $0<c_{0}<c_{1}<\infty$ and $d\geq 1$ such that

We also need the following result regarding the metric entropy of the Grassmannian manifold $G(k,r)$ due to Szarek Szarek82 .

where $c_{0},c_{1}$ are absolute constants. Moreover, for any $\mathbf{V}\in O(k,r)$ and any $\alpha\in(0,1)$,

Proof of Theorem 8 For the purpose of lower bound, we consider the special case of $\lambda_{1}=\cdots=\lambda_{r}=\lambda$, that is, $\bolds{\Sigma}=\lambda\mathbf{V}\mathbf{V}^{\prime}+\mathbf{I}_{k}$. Note that the Kullback–Leibler divergence between normal distributions is given by $D(N(0,\bolds{\Sigma}_{1})||N(0,\bolds{\Sigma}_{0}))=\frac{1}{2}(\operatorname{Tr}(\bolds{\Sigma}_{0}^{-1}\bolds{\Sigma}_{1}-\mathbf{I}_{k})-\log\det\bolds{\Sigma}_{0}^{-1}\bolds{\Sigma}_{1})$. Then for any $\mathbf{U},\mathbf{V}\in O(k,r)$, we have

where the first and second inequalities are by the matrix inversion lemma and the fact that $\operatorname{Tr}(\mathbf{V}\mathbf{V}^{\prime})=\operatorname{Tr}(\mathbf{V}^{\prime}\mathbf{V})=r$, respectively. In view of (47), we have $A=nh(\lambda)/2$. Applying Proposition 3 with $\epsilon_{0}=\sqrt{2(r\wedge(k-r))}$ yields the desired (46).

Proof of Theorem 3 Let $\Theta=\Theta_{0}(s,p,r,\lambda)$. By definition (13), $k^{*}_{0}$ coincides with $s$. In view of the fact that $(a\wedge b)+(c\wedge d)\geq(a\wedge c)(b+d)$, it is sufficient to prove the following inequalities separately:

Inequality (53) follows from an oracle argument: consider the following sub-collection:

Split the data matrix according to $\mathbf{X}=[\mathbf{X}_{1},\mathbf{X}_{2}]$, where $\mathbf{X}_{1}$ consists of the first $s$ columns. Let $\bolds{\Lambda}=\operatorname{diag}({\lambda_{1},\ldots,\lambda_{s}})$. Then the rows of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ are i.i.d. according to $\mathcal{N}(0,\mathbf{V}_{1}\bolds{\Lambda}\mathbf{V}_{1}^{\prime}+\mathbf{I}_{s})$ and $N(0,\mathbf{I}_{p-s})$, respectively. Therefore a sufficient statistic for estimating $\mathbf{V}$ is $\mathbf{X}_{1}$. This reduces the problem to an $s$-dimensional unconstrained PCA problem. Applying the lower bound in Theorem 8 yields (53).

Inequality (54) follows from existing results on rank-one estimation (e.g., Birnbaum12 , Vu12 ). To make the argument rigorous, we focus on the special case where $\{{\mathbf{v}_{2},\ldots,\mathbf{v}_{r}}\}$ are fixed to be standard basis. Denote the following sub-collection:

2 Proof of Theorem 4

We first state a few technical lemmas (proved in Section 7.10 in the supplementary material appsm ) and an oracle upper bound (proved in Section 7.9 in the supplementary material appsm ), which, in view of the lower bound in Theorem 8, gives the optimal rates of the regular PCA problem. Some of the proofs are relegated to the supplementary material appsm .

Let $a,b,c>0$. Then $ax^{2}\leq bx+c$ implies that $x^{2}\leq\frac{b^{2}}{a^{2}}+\frac{2c}{a}$.

Since $|x-\frac{b}{2a}|\leq\frac{\sqrt{b^{2}+4ac}}{2a}$, we have $x^{2}\leq\frac{b^{2}+b^{2}+4ac}{2a^{2}}$.

Let $\bolds{\Sigma}=\mathbf{I}_{p}+\mathbf{V}\mathbf{D}\mathbf{V}^{\prime}$. For any $\mathbf{T}\in O(p,r)$, we have

Let ${X_{1},\ldots,X_{N}}$ be i.i.d. such that

Let $\mathbf{E}$ be a symmetric positive definite matrix. Let $\mathbf{F}$ be a symmetric matrix. Then $|\langle\mathbf{E},\mathbf{F}\rangle|\leq\|{\mathbf{F}}\|\operatorname{Tr}(\mathbf{E})$.

This is a special case of von Neumann’s trace inequality.

Let $\bolds{\Theta}\in\mathcal{F}_{q}(s,p)$ and $k\in[p]$, where $\mathcal{F}_{q}(s,p)$ is defined in (43). Let $\|\bolds{\Theta}_{(i)*}\|$ denote its $i$th largest row norm. Then

By the definition of $\mathcal{F}_{q}(s,p)$ in (43), we have

Let $p=k$ and $r\in[k]$. Let $n\geq C_{0}(r+\log\lambda)$ and $\lambda\geq C_{0}\sqrt{{\log{(n)}}/{n}}$ for some sufficiently large constant $C_{0}$. Let $\widehat{\mathbf{V}}\in O(k,r)$ be formed by the $r$ leading singular vectors of the sample covariance matrix $\mathbf{S}$. Let $\Theta=\Theta_{0}(k,k,r,\lambda,\kappa)$. Then

Proof of Theorem 4 Before delving into the details, we give an outline of the proof as follows: {longlist}[(3)]

We decompose the risk into a summation of three terms, namely the approximation error, oracle risk and excess risk, the first two of which are upper bounded in Lemma 7 and Theorem 9, respectively.

The excess risk is controlled by a careful concentration-of-measure analysis, which forms the core of the proof. We also remark that by (8), (13) and condition (25), we have

To see this, first note that $k_{0}^{*}\geq r$ by (8) directly. When $q\in(0,2)$, if $k_{q}^{*}=p$, then $k_{q}^{*}\geq r$. Otherwise, we have

Here the first inequality comes from (13), the second is due to condition (25), the third holds since $k_{q}^{*}\geq 1$ and the last holds for sufficiently large $C_{0}$ in view of (8). {longlist}

. Fix $\mathbf{V}\in O(p,r)\cap\mathcal{F}_{q}(s,p)$. We assume that $q>0$. Note that this step is superfluous if $q=0$ since $\mathbf{V}$ is already sparse. Let $k=k_{q}^{*}$ be defined in (13). Let $\mathcal{B}(k)=\{B\subset[p]\dvtx|B|=k\}$. Let $A\in\mathcal{B}(k)$ denote the collection of row indices of $\mathbf{V}$ corresponding to the $k$ largest row norm. Put

Put $\mathbf{U}=\mathbf{J}_{A}\mathbf{V}$. Then

where (65) follows from applying Lemma 7, (66) follows from the choice of $k=k^{*}_{q}$ in (13), and (67) is implied by the assumption (25). Therefore

By definition of the maximizer $B^{*}$ in (22), $\langle\mathbf{S}_{(2)},\mathbf{V}_{A}\mathbf{V}_{A}^{\prime}-\mathbf{V}_{*}\mathbf{V}_{*}^{\prime}\rangle\leq 0$. In view of Lemma 3, we have

The hard part is to control the third term (the worst-case fluctuation) in (71). To this end, we decompose the sample covariance matrix as

We first deal the inner product with $\mathbf{G}$: write $\langle\mathbf{G},\widehat{\mathbf{V}}_{A}\widehat{\mathbf{V}}_{A}^{\prime}-\widehat{\mathbf{V}}_{*}\widehat{\mathbf{V}}_{*}^{\prime}\rangle=\langle\mathbf{G},\widehat{\mathbf{V}}_{A}\widehat{\mathbf{V}}_{A}^{\prime}-\mathbf{V}\mathbf{V}^{\prime}\rangle-\langle\mathbf{G},\widehat{\mathbf{V}}_{*}\widehat{\mathbf{V}}_{*}^{\prime}-\mathbf{V}\mathbf{V}^{\prime}\rangle$. Note that

where (76) is due to (19) and (75) is a consequence of Lemma 6, in view of the fact that $\mathbf{I}_{r}-\mathbf{V}^{\prime}\widehat{\mathbf{V}}_{A}\widehat{\mathbf{V}}_{A}^{\prime}\mathbf{V}$ is symmetric positive semi-definite while $\mathbf{D}(\frac{1}{n}\mathbf{U}_{(2)}^{\prime}\mathbf{U}_{(2)}-\mathbf{I}_{r})\mathbf{D}$ is symmetric. Similarly, we have

which has zero trace and unit Frobenius norm. Recall that $\widehat{\mathbf{V}}_{*}=\widehat{\mathbf{V}}_{B^{*}}$. Then

Assembling (72), (78) and (6.2), we can upper bound the excess risk by

Now we combine the risk decomposition (71) with the upper bounds above to control the risk of our aggregated estimator $\widehat{\mathbf{V}}_{*}$: to simplify notation, denote

Introduce the event $E=\{M\leq\frac{1}{24\kappa}\}$. By assumption (26), $r\leq c^{\prime\prime}n$ for a sufficiently small constant $c^{\prime\prime}$. Then there exists a constant $c^{\prime}>0$ only depending on $\kappa$, such that $\frac{1}{24\kappa}\geq 2(\sqrt{r\over n}+t)+(\sqrt{r\over n}+t)^{2}$, where $t=\sqrt{\frac{\log(c^{\prime}nh(\lambda))}{n}}$. Applying Proposition 4 in the supplementary material appsm yields

Conditioning on the event $E$ and using Lemma 2, we have

Recall from (19) that the loss function is upper bounded by $r\wedge(p-r)$. Taking expectation on both sides of (84), and using (83) together with the Cauchy–Schwarz inequality, we have

In view of the oracle upper bound in Theorem 9, we have

By (69), if $q>0$, the approximation is upper bounded by

If $q=0$, then $\Delta=0$. To control the right-hand side of (86), it boils down to upper bound $\mathsf{E}T^{2}$. In the sequel we shall prove that

for some absolutely constant $C$. Plugging (87), (88) and (89) into (86), we arrive at

where the constant $C^{\prime}$ only depends on $\kappa$. In the special case of $q=0$, the approximation error is $\Delta=0$, which implies that the second term in (91) is zero. Hence we have the following stronger result:

where $\Psi_{0}$ is defined in (12). Then (91) and (6.2) imply the statement of the theorem for $q>0$ and $q=0$, respectively.

To finish the proof of the theorem, it remains to establish (89). To this end, recall that $\mathbf{K}_{B}$ is symmetric and $\operatorname{Tr}(\mathbf{K}_{B})=0$. By the definitions of $T$ and $\mathbf{H}$ in (6.2) and (74), respectively, we have

Assembling (93) with (96)–(95) and using the fact that $(a+b)^{2}\leq 2(a^{2}+b^{2})$, we arrive at

where we used $\frac{k}{n}\log\frac{p}{k}\leq 1$ implied by the assumption (26).

It then remains to establish (96)–(97). Note that the collection $\{\mathbf{K}_{B}\dvtx B\in\mathcal{B}(k)\}$ belongs to the $\sigma$-algebra generated by the first sample $\mathbf{X}_{(1)}$, which is independent of $(\mathbf{Z}_{(2)},\mathbf{U}_{(2)})$. By conditioning on $\mathbf{X}_{(1)}$, we can treat $\{\mathbf{K}_{B}\dvtx B\in\mathcal{B}(k)\}$ as fixed matrices. ∎\noqed

Proof of (96) For each fixed $B\in\mathcal{B}(k)$, $\mathbf{K}_{B}\perp\!\!\!\perp\mathbf{Z}_{(2)}$. Applying Lemma 4, we have

for some $W\sim N(0,1)$ independent of $\mathbf{U}_{(2)}$.

Consequently, $\langle\mathbf{V}\mathbf{D}\mathbf{U}_{(2)}^{\prime}\mathbf{Z}_{(2)},\mathbf{K}_{B}\rangle$ is stochastically dominated by$\sqrt{\lambda_{1}}\|{\mathbf{U}_{(2)}}\||W|$. Since $\mathbf{U}_{(2)}$ is an $n\times r$ standard Gaussian matrix, Lemma 10 in the supplementary material appsm yields

which the last inequality follows from (102) and the Chernoff bound $\mathsf{P}(W\geq\sqrt{2}t)\leq\frac{1}{2}\exp(-t)$. Therefore,

Applying Lemma 5 with $N={p\choose k}$ yields

which, in view of $r\leq k$, implies the desired (97).

3 Proof of Theorem 7

We prove the theorem in three steps. First, we verify that the “whitening” procedure in step 3 of the reduction scheme can be performed. Next, we investigate the signal-to-noise ratio of the regression problem conditional on the values of $\mathbf{U}$ and $\mathbf{Z}^{0}$. Finally, we derive the desired rates by using Theorem 6 and Wedin’s sin-theta theorem Wedin72 . {longlist}[($2^{\circ}$)]

As a first step, we verify that the “whitening” step is indeed possible, which requires that $\sigma_{r}(\mathbf{B})>0$. To this end, let $J=\operatorname{supp}(\mathbf{V}^{0})$. Since $\mathbf{B}=\mathbf{U}\mathbf{D}\mathbf{V}^{\prime}\mathbf{V}^{0}+\mathbf{Z}^{0}\mathbf{V}^{0}$, we have

By our assumption on $\mathbf{V}^{0}$, condition (41) is satisfied with probability at least $1-C/[nh(\lambda)]$. By Lemma 10 in the supplementary material appsm and the union bound,

holds with probability at least $1-C/[nh(\lambda)]$. Note that assumption (26) implies that $n\geq C_{0}r$ and that $n\geq C_{0}\log[nh(\lambda)]$. Thus, for sufficiently large $C_{0}$ in (26), the first inequality in (6.3) leads to $\sigma_{r}(\mathbf{U})\geq\frac{2}{3}\sqrt{n}$. Together with $\sigma_{r}(\mathbf{D})=\sqrt{\lambda_{r}}$, the first term in (6.3) is thus lower bounded by $\frac{1}{3}\sqrt{n\lambda_{r}}$, and hence

with probability at least $1-C/[nh(\lambda)]$.

Turning to the second term in (6.3), we first note that it is upper bounded by $\max_{I\subset[p],|I|=k^{*}_{q}}\|\mathbf{Z}^{0}_{I}\|$ conditioned on the event that $|J|\leq k_{q}^{*}$. Note that for any $t>0$, we have

with probability at least $1-C/[nh(\lambda)]$, where the last inequality holds because the assumption (26) implies that $k_{q}^{*}\leq n/4$ and $t^{*}\leq n/2$ as long as $C_{0}$ is sufficiently large.

Under the assumption that $\lambda_{r}\geq C_{0}$ for some sufficiently large $C_{0}>36$, (105) and (106) lead to $\sigma_{r}(\mathbf{B})\geq c\sqrt{n\lambda_{r}}>0$ with probability at least $1-C/[nh(\lambda)]$. This completes the first step in the proof.

Let ${\bar{\mathbf{A}}}=\tfrac{1}{\sqrt{2}}\mathbf{A}\mathbf{R}\mathbf{C}^{-1}=\tfrac{1}{\sqrt{2}}\mathbf{D}\mathbf{U}^{\prime}\mathbf{B}\mathbf{R}\mathbf{C}^{-1}=\tfrac{1}{\sqrt{2}}\mathbf{D}\mathbf{U}^{\prime}\mathbf{L}$. Then $\bolds{\Theta}=\mathbf{V}{\bar{\mathbf{A}}}$ in (32). In the second step, we show that there exist two constants $C_{2}>C_{1}>0$ depending only on $\kappa$, such that with probability at least $1-C/[nh(\lambda)]$,

To this end, note that (6.3) and assumption (26) imply

holds with probability at least $1-C/[nh(\lambda)]$. Under the same assumption, Lemma 10 in the supplementary material appsm implies

Next we show that, conditioned on the event that (107) holds, the signal matrix $\bolds{\Theta}$ lies in $\mathcal{F}_{q}(s^{\prime},p)$ where

where the middle inequality is due to (107), the last inequality follows from the assumption that $\lambda\geq C_{0}$ and the first inequality is due to $\|\bolds{\Theta}\|_{q,w}\leq\|\mathbf{V}\|_{q,w}\|\bar{\mathbf{A}}\|^{q}$, which is a consequence of equation (110) in Section 7.1 of the supplementary material appsm .

Let $k^{\prime}$ be defined in (44). We show that whenever (108) holds, we have

Let $E$ denote the event that both (6.3) and (107) hold. Then

Here, the last inequality holds because the loss function is upper bounded by $r$ and $\mathsf{P}(E^{c})\leq C/[nh(\lambda)]$.

To further bound the first term on the rightmost hand side, we note that $E$ is completely determined by $\mathbf{U}$ and $\mathbf{Z}^{0}$. Hence, it is nonrandom conditioned on $\mathbf{U}$ and $\mathbf{Z}^{0}$. Thus

Supplement to “Sparse PCA: Optimal rates and adaptive estimation” \slink[doi]10.1214/13-AOS1178SUPP \sdatatype.pdf \sfilenameaos1178_supp.pdf \sdescriptionWe provide proofs for all the remaining theoretical results in the paper. The proofs rely on results in cover , Davidson01 , Davis70 , Johnstone01ss , KT59 , Laurent00 and Tsybakov09 .

References