Sparse principal component analysis and iterative thresholding

Introduction

In many contemporary datasets, if we organize the $p$-dimensional observations $x_{1},\ldots,x_{n}$, into the rows of an $n\times p$ data matrix $X$, the number of features $p$ is often comparable to, or even much larger than, the sample size $n$. For example, in biomedical studies, we usually have measurements on the expression levels of tens of thousands of genes, but only for tens or hundreds of individuals. One of the crucial issues in the analysis of such “large $p$” datasets is dimension reduction of the feature space.

As a classical method, principal component analysis (PCA) pear01 , hote33 reduces dimensionality by projecting the data onto the principal subspace spanned by the $m$ leading eigenvectors of the population covariance matrix $\Sigma$, which represent the principal modes of variation. In principle, one expects that for some $m<p$, most of the variance in the data is captured by these $m$ modes. Thus, PCA reduces the dimensionality of the feature space while retaining most of the information in data. In addition, projection to a low-dimensional space enables visualization of the data. In practice, $\Sigma$ is unknown. Classical PCA then estimates the leading population eigenvectors by those of the sample covariance matrix $S$. It performs well in the traditional data setting where $p$ is small and $n$ is large ande63 .

In high-dimensional settings, a collection of data can be modeled by a low-rank signal plus noise structure, and PCA can be used to recover the low-rank signal. In particular, each observation vector $x_{i}$ can be viewed as an independent instantiation of the following generative model:

Here, $\mu$ is the mean vector, $A$ is a $p\times\bar{m}$ deterministic matrix of factor loadings, $u_{i}$ is an $\bar{m}$-vector of random factors, $\sigma>0$ is the noise level and $z_{i}$ is a $p$-vector of white noise. For instance, in chemometrics, $x_{i}$ can be a vector of the logarithm of the absorbance or reflectance spectra measured with noise, where the columns of $A$ are characteristic spectral responses of different chemical components, and $u_{i}$’s the concentration levels of these components vafi09 . The number of observations are relatively few compared with the number of frequencies at which the spectra are measured. In econometrics, $x_{i}$ can be the returns for a collection of assets, where the $u_{i}$’s are the unobservable random factors tsay05 . The assumption of additive white noise is reasonable for asset returns with low frequencies (e.g., monthly returns of stocks). Here, people usually look at tens or hundreds of assets simultaneously, while the number of observations are also at the scale of tens or hundreds. In addition, model (1) represents a big class of signal processing problems wk85 . Without loss of generality, we assume $\mu=0$ from now on.

In this paper, our primary interest lies in PCA of high-dimensional data generated as in (1). Let the covariance matrix of $u_{i}$ be $\Phi$ which is of full rank. Suppose that $A$ has full column rank and that $u_{i}$ and $z_{i}$ are independent. Then the covariance matrix of $x_{i}$ becomes

Here, $\lambda_{1}^{2}\geq\cdots\geq\lambda_{\bar{m}}^{2}>0$ are the eigenvalues of $A\Phi A^{\prime}$, with $q_{j}$, $j=1,\ldots,\bar{m}$, the associated eigenvectors. Therefore, the $j$th eigenvalue of $\Sigma$ is $\lambda_{j}^{2}+\sigma^{2}$ for $j=1,\ldots,\bar{m}$, and $\sigma^{2}$ otherwise. Since there are $\bar{m}$ spikes $(\lambda_{1}^{2},\ldots,\lambda_{\bar{m}}^{2})$ in the spectrum of $\Sigma$, (2) has been called the spiked covariance model in the literature john01 . Note that we use $\lambda_{j}^{2}$ to denote the spikes rather than $\lambda_{j}$ used previously in the literature pajo07 . For data with such a covariance structure, it makes sense to project the data onto the low-dimensional subspaces spanned by the first few $q_{j}$’s. Here and after, $\bar{m}$ denotes the number of spikes in the model, and $m$ is the target dimension of the principal subspace to be estimated, which is no greater than $\bar{m}$.

Classical PCA encounters both practical and theoretical difficulties in high dimensions. On the practical side, the eigenvectors found by classical PCA involve all the $p$ features, which makes their interpretation challenging. On the theoretical side, the sample eigenvectors are no longer always consistent estimators. Sometimes, they can even be nearly orthogonal to the target direction. When both $n,p\to\infty$ with $n/p\to c\in(0,\infty)$, at different levels of rigor and generality, this phenomenon has been examined by a number of authors reva96 , lu02 , hora04 , paul07 , nadl08 , onat09 under model (2). See juma09 for similar results when $p\to\infty$ and $n$ is fixed.

In recent years, to facilitate interpretation, researchers have started to develop sparse PCA methodologies, where they seek a set of sparse vectors spanning the low-dimensional subspace that explains most of the variance. See, for example, jotrud03 , zohati06 , dagh07 , shhu08 , solo08 , witiha09 . These approaches typically start with a certain optimization formulation of PCA and then induce a sparse solution by introducing appropriate penalties or constraints.

On the other hand, when $\Sigma$ indeed has sparse leading eigenvectors in the current basis (perhaps after transforming the data), it becomes possible to estimate them consistently under high-dimensional settings via new estimation schemes. For example, under normality assumption, when $\Sigma$ only has a single spike, that is, when $\bar{m}=1$ in (2), Johnstone and Lu jolu09 proved consistency of PCA obtained on a subset of features with large sample variances when the leading eigenvalue is fixed and $(\log{p})/n\to 0$. Under the same single spike model, if in addition the leading eigenvector has exactly $k$ nonzero loadings, Amini and Wainwright amwa09 studied conditions for recovering the nonzero locations using the methods in jolu09 and dagh07 , and Shen et al. shshma11 established conditions for consistency of a sparse PCA method in shhu08 when $p\to\infty$ and $n$ is fixed. For the more general multiple component case, Paul and Johnstone pajo07 proposed an augmented sparse PCA method for estimating each of the leading eigenvectors, and showed that their procedure attains near optimal rate of convergence under a range of high-dimensional sparse settings when the leading eigenvalues are comparable and well separated. Notably, these methods all focus on estimating individual eigenvectors.

In this paper, we focus primarily on finding principal subspaces of $\Sigma$ spanned by sparse leading eigenvectors, as opposed to finding each sparse vector individually. One of the reasons is that individual eigenvectors are not identifiable when some leading eigenvalues are identical or close to each other. Moreover, if we view PCA as a dimension reduction technique, it is the low-dimensional subspace onto which we project data that is of the greatest interest.

The contribution of the current paper is threefold. First, we propose to estimate principal subspaces. This is natural for the purpose of dimension reduction and visualization, and avoids the identifiability issue for individual eigenvectors. Second, we construct a new algorithm to estimate the subspaces, which is efficient in computation and easy to implement. Last but not least, we derive convergence rates of the resulting estimator under the spiked covariance model when the eigenvectors are sparse.

The rest of the paper is organized as follows. In Section 2, we frame the principal subspace estimation problem and propose the iterative thresholding algorithm. The statistical properties and computational complexity of the algorithm are examined in Sections 3 and 4 under normality assumption. Simulation results in Section 5 demonstrate its competitive performance. Section 6 presents the proof of the main theorems.

Reproducible code: The Matlab package SPCALab implementing the proposed method and producing the tables and figures of the current paper is available at the author’s website.

Methodology

We use $C,C_{0},C_{1}$, etc. to represent constants, though their values might differ at different occurrences. For real numbers $a$ and $b$, let $a\vee b=\max(a,b)$ and $a\wedge b=\min(a,b)$. We write $a_{n}=O(b_{n})$, if there is a constant $C$, such that $|a_{n}|\leq Cb_{n}$ for all $n$, and $a_{n}=o(b_{n})$ if $a_{n}/b_{n}\to 0$ as $n\to\infty$. Moreover, we write $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. Throughout the paper, we use $\nu$ as the generic index for features, $i$ for observations, $j$ for eigenvalues and eigenvectors and $k$ for iterations in the algorithm to be proposed.

2 Framing the problem: Principal subspace estimation

To measure the accuracy of an estimator $\widehat{\mathcal{S}}$ for a subspace $\mathcal{S}$, note that each linear subspace is associated with a unique projection matrix onto it. Let $P$ and $\widehat{P}$ be the projection matrices associated with $\mathcal{S}$ and $\widehat{\mathcal{S}}$, respectively. The distance between $\mathcal{S}$ and $\widehat{\mathcal{S}}$ is given by the spectral norm of the difference between $P$ and $\widehat{P}$: $\operatorname{dist}(\mathcal{S},\widehat{\mathcal{S}})=\|P-\widehat{P}\|$; see gova96 , Section 2.6.3. Thus, we can define a loss function by the squared distances between $\mathcal{S}$ and $\widehat{\mathcal{S}}$,

By definition, this loss function measures the maximum possible discrepancy between the projections of any unit vector onto the two subspaces. The loss ranges in $$, and equals zero if and only if$\widehat{\mathcal{S}}=\mathcal{S}$. When$\operatorname{dim}(\widehat{\mathcal{S}})\neq\operatorname{dim}(\mathcal{S})$, we have$L(\mathcal{S},\widehat{\mathcal{S}})=1$. Geometrically, it equals the squared sine of the largest canonical angle between$\mathcal{S}$and$\widehat{\mathcal{S}}$ (stsu90 , Theorem 5.5). Throughout the paper, we use the loss function (3) for principal subspace estimation.

3 Orthogonal iteration

Given a positive definite matrix $A$, a standard technique to compute its leading eigenspace is orthogonal iteration gova96 . When only the first eigenvector is sought, it is also known as the power method.

To state the orthogonal iteration method, we note that for any $p\times m$ matrix $T$, when $p\geq m$, we could decompose it into the product of two matrices $T=QR$, where $Q$ is $p\times m$ orthonormal and $R$ is $m\times m$ upper triangular. This decomposition is called QR factorization and can be computed using Gram–Schmidt orthogonalization and other numerical methods gova96 . Suppose $A$ is $p\times p$, and we want to compute its leading eigenspace of dimension $m$. Starting with a $p\times m$ orthonormal matrix $Q^{(0)}$, orthogonal iteration generates a sequence of $p\times m$ orthonormal matrices $Q^{(k)}$, $k=1,2,\ldots\,$, by alternating the following two steps till convergence: {longlist}[(2)]

QR factorization: $Q^{(k)}R^{(k)}=T^{(k)}$. Denote the orthonormal matrix at convergence by $Q^{(\infty)}$. Then its columns are the leading eigenvectors of $A$, and $\operatorname{ran}(Q^{(\infty)})$ gives the eigenspace. In practice, one terminates the iteration once $\operatorname{ran}(Q^{(k)})$ stabilizes.

When we apply orthogonal iteration directly to the sample covariance matrix $S$, it gives the classical PCA result, which could be problematic in high dimensions. Observe that all the $p$ features are included in orthogonal iteration. When the dimensionality is high, not only the interpretation is hard, but the variance accumulated across all the features becomes so high that it makes consistent estimation impossible.

If the eigenvectors spanning $\mathcal{P}_{m}$ are sparse in the current basis, one sensible way to reduce estimation error is to focus only on those features at which the leading eigenvectors have large values, and to estimate other features by zeros. Of course, one introduces bias this way, but hopefully it is much smaller compared to the amount of variance thus reduced.

The above heuristics lead to the estimation scheme in the next subsection which incorporates this feature screening idea in orthogonal iteration.

4 Iterative thresholding algorithm

Let $S=\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\prime}$ be the sample covariance matrix. An effective way to incorporate feature screening into orthogonal iteration is to “kill” small coordinates of the $T^{(k)}$ matrix after each multiplication step, which leads to the estimation scheme summarized in Algorithm 1. Although the later theoretical study is conducted under normality assumption, Algorithm 1 itself is not confined to normal data.

In addition to the two basic orthogonal iteration steps, Algorithm 1 adds a thresholding step in between them, where we threshold each element of $T^{(k)}$ with a user-specified thresholding function $\eta$ which satisfies

The ranges of $\widehat{Q}{}^{(k)}$ and $\widehat{T}{}^{(k)}$ are the same because QR factorization only amounts to a basis change within the same subspace. However, as in orthogonal iteration, the QR step is essential for numerical stability, and should not be omitted. Moreover, although the algorithm is designed for subspace estimation, the column vectors of $\widehat{Q}{}^{(\infty)}$ can be used as estimators of leading eigenvectors.

Algorithm 1 requires an initial orthonormal matrix $\widehat{Q}{}^{(0)}$. It can be generated from the “diagonal thresholding” sparse PCA algorithm jolu09 . Its multiple eigenvector version is summarized in Algorithm 2. Here, for any set $I$, $\operatorname{card}(I)$ denotes its cardinality.

Given the output $\widehat{Q}_{B}=[\widehat{q}_{1},\ldots,\widehat{q}_{\operatorname{card}(B)}]$ of Algorithm 2, we take $\widehat{Q}{}^{(0)}=[\widehat{q}_{1},\ldots,\widehat{q}_{m}]$. When $\sigma^{2}$ is unknown, we could replace it by an estimator $\widehat{\sigma}^{2}$ in the definition of $B$. For example, for normal data, Johnstone and Lu jolu09 suggested

When available, subject knowledge could also be incorporated into the construction of $\widehat{Q}{}^{(0)}$. Algorithm 1 also requires inputs for the $\gamma_{nj}$’s and subspace dimension $m$. Under normality assumption, we give explicit specification for them in (8) and (20) later. Under the conditions of the later Section 3, $B$ is nonempty with probability tending to $1$, and so $\widehat{Q}{}^{(0)}$ is well defined.

Convergence

For normal data, to obtain the error rates in later Theorems 3.1 and 3.2, we can terminate Algorithm 1 after $K_{s}$ iterations with $K_{s}$ given in (9). In practice, one could also stop iterating if the difference between successive iterates becomes sufficiently small, for example, when $L(\operatorname{ran}(\widehat{Q}{}^{(k)}),\operatorname{ran}(\widehat{Q}{}^{(k+1)}))\leq n^{-2}$. We suggest this empirical stopping rule because $n^{-2}$ typically tends to zero faster than the rates we shall obtain, and so intuitively it should not change the statistical performance of the resulting estimator. In simulation studies reported in Section 5, the difference in numerical performance between the outputs based on this empirical stopping rule and those based on the theoretical rule (9) is negligible compared to the estimation errors. Whether Algorithm 1 always converges numerically is an interesting question left for possible future research.

Bibliographical note

When $m=1$, Algorithm 1 is similar to the algorithms proposed in shhu08 , witiha09 and yuzh11 . When $m>1$, all these methods propose to iteratively find the first leading eigenvectors of residual covariance matrices, which becomes different from our approach.

Statistical properties

This section is devoted to analyzing the statistical properties of Algorithm 1 under normality assumption. After some preliminaries, we first establish the convergence rates for subspace estimation in a special yet interesting case in Section 3.1. Then we introduce a set of general assumptions in Section 3.2 and a few key quantities in Section 3.3. Section 3.4 states the main results, which include convergence rates for principal subspace estimation under general assumptions and a correct exclusion property. In addition, we derive rates for estimating individual eigenvectors. For conciseness, we first state all the results assuming a suitable target subspace dimension $m\leq\bar{m}$ is given. In Section 3.5, we discuss how to choose $m$ and estimate $\bar{m}$ based on data.

We start with some preliminaries. Under normality assumption, $x_{1},\ldots,x_{n}$ are i.i.d. $N_{p}(0,\Sigma)$ distributed, with $\Sigma$ following model (2). Further assume $\sigma^{2}$ is known—though this assumption could be removed by estimating $\sigma^{2}$ using, say, $\widehat{\sigma}^{2}$ in (5). Since one can always scale the data first, we assume $\sigma^{2}=1$ from now on. Thus, (1) reduces to the orthogonal factor form

Here, $v_{ij}$ are i.i.d. standard normal random factors, which are independent of the i.i.d. white noise vectors $z_{i}\sim N_{p}(0,I)$, and $\{q_{j},1\leq j\leq\bar{m}\}$ is a set of leading eigenvectors of $\Sigma$. In what follows, we use $n$ to index the size of the problem. So the dimension $p=p(n)$ and the spikes $\lambda_{j}^{2}=\lambda_{j}^{2}(n)$ can be regarded as functions of $n$, while both $\bar{m}$ and $m$ remain fixed as $n$ grows.

Let $p_{n}=p\vee n$. We obtain the initial matrix $\widehat{Q}{}^{(0)}$ in Algorithm 1 by applying Algorithm 2 with

In Algorithm 1, the threshold levels are set at

To facilitate understanding, we first state the convergence rates for principal subspace estimation in a special case.

Recall that $h(x)=x^{2}/(x+1)$. Under the above setup, we have the following upper bound for subspace estimation error.

Under the above setup, for sufficiently large constants $\alpha,\gamma{>2\sqrt{3}}$ in (7) and (8), there exist constants $C_{0}$, $C_{1}=C_{1}(\gamma,r,m)$ and $C_{2}$, such that for sufficiently large $n$, uniformly over all $\Sigma$ with $\|q_{j}\|_{r}\leq s$ for $1\leq j\leq\bar{m}$, with probability at least $1-C_{0}p_{n}^{-2}$, we have $K_{s}\asymp\log{n}$ and the subspace estimator $\widehat{\mathcal{P}}_{m}^{({K_{s}})}=\operatorname{ran}(\widehat{Q}{}^{({K_{s}})})$ of Algorithm 1 satisfies

where $g_{m}(\lambda)=\frac{(\lambda_{1}^{2}+1)(\lambda_{m+1}^{2}+1)}{(\lambda_{m}^{2}-\lambda_{m+1}^{2})^{2}}$.

The upper bound in Theorem 3.1 consists of two terms. The first is a “nonparametric” term, which can be decomposed as the product of two components. The first component, $\bar{m}s^{r}[nh(\lambda_{m}^{2})/\log{p}]^{r/2}$, up to a multiplicative constant, bounds the number of coordinates used in estimating the subspace, while the second component, $\log{p}/[nh(\lambda_{m}^{2})]$, gives the average error per coordinate. The second term in the upper bound, $g_{m}(\lambda)(\log{p})/n$, up to a logarithmic factor, has the same form as the cross-variance term in the “fixed $p$, large $n$” asymptotic limit for classical PCA; cf. ande63 , Theorem 1. We call it a “parametric” error term, because it always arises when we try to separate the first $m$ eigenvectors from the rest, regardless of how sparse they are. Under the current setup, both terms converge to as $n\to\infty$, which establishes the consistency of our estimator.

Since $\alpha_{n}$ and $\gamma_{nj}$ and the stopping rule (9) do not involve any unknown parameter, the theorem establishes the adaptivity of our estimator: the optimal rate of convergence in Theorem 3.1 is obtained without any knowledge of the power $r$, the radius $s$ or the spikes $\lambda_{j}^{2}$. Last but not least, the estimator could be obtained in $O(\log{n})$ iterations and holds for all thresholding function $\eta$ satisfying (4).

Later in Section 3.4, Theorem 3.2 establishes analogous convergence rates, but for a much wider range of high-dimensional sparse settings. In particular, the above result will be extended simultaneously along two directions: {longlist}[(2)]

the spikes $\lambda_{1}^{2},\ldots,\lambda_{\bar{m}}^{2}$ will be allowed to scale as $n\to\infty$, and $\lambda_{m+1}^{2},\ldots,\break\lambda_{\bar{m}}^{2}$ could even be of smaller order as compared to the first $m$ spikes;

2 Assumptions

We now state assumptions for the general theoretical results in Section 3.4.

As outlined above, the first extension of the special case is to allow the spikes $\lambda_{j}^{2}=\lambda_{j}^{2}(n)>0$ to change with $n$, though the dependence will usually not be shown explicitly. Recall that $p_{n}=p\vee n$; we impose the following growth rate condition on $p$ and the $\lambda_{j}^{2}$’s.

As $n\to\infty$, we have: {longlist}[(3)]

the dimension $p$ satisfies $(\log{p})/n=o(1)$;

the largest spike $\lambda_{1}^{2}$ satisfies $\lambda_{1}^{2}=O(p_{n})$; the smallest spike $\lambda_{\bar{m}}^{2}$ satisfies $\log(p_{n})=o(n\lambda_{\bar{m}}^{4})$; and their ratio satisfies ${\lambda_{1}^{2}}/{\lambda_{\bar{m}}^{2}}=O(n[\log(p_{n})/n]^{{1/2}+r/4})$;

$\lim_{n\to\infty}\lambda_{1}^{2}/(\lambda_{j}^{2}-\lambda_{j+1}^{2})\in[1,\infty]$ exists for $j=1,\ldots,\bar{m}$, with $\lambda^{2}_{\bar{m}+1}=0$.

The first part of Condition GR requires the dimension to grow at a sub-exponential rate of the sample size. The second part ensures that the spikes grow at most at linear rate with $p_{n}$, and are all of larger magnitude than $\sqrt{\log(p_{n})/n}$. In addition, the condition on the ratio $\lambda_{1}^{2}/\lambda_{\bar{m}}^{2}$ allows us to deal with the interesting cases where the first several spikes scale at a faster rate with $n$ than the others. This is more flexible than the assumption previously made in pajo07 that all the spikes grow at the same rate. The third part requires $\lim_{n\to\infty}\lambda_{1}^{2}/(\lambda_{j}^{2}-\lambda_{j+1}^{2})$ to exist for each $1\leq j\leq\bar{m}$, but the limit can be infinity.

For general results, we allow the radii $s_{j}$’s to depend on or even diverge with $n$, though we require that they do not grow too rapidly, so the leading eigenvectors are indeed sparse. This leads to the following sparsity condition.

This type of condition also appeared in a previous study of individual eigenvector estimation in the multiple component spiked covariance model pajo07 . The condition is, for example, satisfied if Condition GR holds and the largest spike $\lambda_{1}^{2}$ is bounded away from zero while the radii $s_{j}$’s are all bounded above by an arbitrarily large constant. That is, if there exists a constant $C>0$, such that $\lambda_{1}^{2}\geq 1/C$ and $s_{j}\leq C$ for all $j\leq\bar{m}$ and all $n$.

It is straightforward to verify that Conditions GR and SP are satisfied by the special case in Section 3.1. We conclude this part with an example.

3 Key quantities

We now introduce a few key quantities which appear later in the general theoretical results.

The first quantity gives the rate at which we distinguish high from low signal coordinates. Recall that $h(x)=x^{2}/(x+1)$. For $j=1,\ldots,\bar{m}$, define

According to paul05 , up to a logarithmic factor, $\tau_{nj}^{2}$ can be interpreted as the average error per coordinate in estimating an eigenvector with eigenvalue $\lambda_{j}^{2}+1$. Thus, a coordinate can be regarded as of high signal if at least one of the leading eigenvectors is of larger magnitude on this coordinate compared to $\tau_{nj}$. Otherwise, we call it a low signal coordinate. We define $H(\beta)$ to be the set of high signal coordinates

Here, $\beta$ is a constant not depending on $n$, the actual value of which will be specified in Theorem 3.2. If $\bar{m}=1$ and $q_{1}$ has $k$ nonzero entries all equal to $1/\sqrt{k}$, then $H$ contains exactly these $k$ coordinates when $k<nh(\lambda_{1}^{2})/[\beta^{2}\log(p_{n})]$, which is guaranteed under Condition SP. In addition, let $L=\{1,\ldots,p\}\setminus H$ be the complement of $H$. Here, $H$ stands for “high,” and $L$ for “low” (also recall $B$ in Algorithm 2, where $B$ stands for “big”). The dependence of $H$, $L$ and $B$ on $n$ is suppressed for notational convenience.

To understand the convergence rate of the subspace estimator stated later in (16), it is important to have an upper bound for $\operatorname{card}(H)$, the cardinality of $H$. To this end, define

The following lemma shows that a constant multiple of $M_{n}$ bounds $\operatorname{card}(H)$. The proof of the lemma is given in supp . Thus, in the general result, $M_{n}$ plays the same role as the term $\bar{m}s^{r}[nh(\lambda^{2})/\log{p}]^{r/2}$ has played in Theorem 3.1.

For sufficiently large $n$, the cardinality of $H=H(\beta)$ satisfies $\bar{m}\leq\operatorname{card}(H)\leq CM_{n}$ for a constant $C$ depending on $\beta$ and $r$.

The last quantity we introduce is related to the “parametric” term in the convergence rate. Let $\lambda^{2}_{\bar{m}+1}=0$. For $j=1,\ldots,\bar{m}$, define

So the second term of the upper bound in Theorem 3.1 is $C_{2}\varepsilon_{nm}^{2}$. For the interpretation of this quantity, we refer to the discussion after Theorem 3.1.

4 Main results

We turn to the statement of main theoretical results.

A key condition for the results is the asymptotic distinguishability (AD) condition introduced below. Recall that all the spikes $\lambda_{j}^{2}$ (hence all the leading eigenvalues) are allowed to depend on $n$. The condition AD will guarantee that the largest few eigenvalues are asymptotically well separated from the rest of the spectrum, and so the corresponding principal subspace is distinguishable.

We say that condition $\operatorname{AD}(j,\kappa)$ is satisfied with constant $\kappa$, if there exists a numeric constant $\kappa\geq 1$, such that for sufficiently large $n$, the gap between the $j$th and the $(j+1)$th eigenvalues satisfies

We define $\operatorname{AD}(0,\kappa)$ and $\operatorname{AD}(\bar{m},\kappa)$ by letting $\lambda_{0}^{2}=\infty$, and $\lambda^{2}_{\bar{m}+1}=0$. So $\operatorname{AD}(0,\kappa)$ holds for any $\kappa\geq 1$. Note that there is always some $1\leq j\leq\bar{m}$ such that condition $\operatorname{AD}(j,\kappa)$ is satisfied. For instance, $\operatorname{AD}(j,\kappa)$ is satisfied with some $\kappa$ for the largest $j$ such that $\lambda_{j}^{2}\asymp\lambda_{1}^{2}$. When the spikes do not change with $n$, condition $\operatorname{AD}(\bar{m},\kappa)$ is satisfied with any constant $\kappa\geq\lambda_{1}^{2}/\lambda_{\bar{m}}^{2}$.

Recall definitions (7)–(9) and (11)–(14). The following theorem establishes the rate of convergence of the principal subspace estimator obtained via Algorithm 1 under relaxed assumptions, which generalizes Theorem 3.1.

Suppose Conditions GR and SP hold, and condition $\operatorname{AD}(m,\kappa)$ is satisfied with some constant $\kappa\geq 1$ for the given subspace dimension $m$. Let the constants $\alpha,\gamma>2\sqrt{3}$ in (7) and (8), and for $c=0.9(\gamma-2\sqrt{3})$, let $\beta=c/\sqrt{m}$ in $H$ (12). Then, there exist constants $C_{0}$, $C_{1}=C_{1}(\gamma,r,m,\kappa)$ and $C_{2}$, such that for sufficiently large $n$, uniformly over $\mathcal{F}_{n}$, with probability at least $1-C_{0}p_{n}^{-2}$, $K_{s}\in[K,2K]$ for

and the subspace estimator $\widehat{\mathcal{P}}_{m}^{({K_{s}})}=\operatorname{ran}(\widehat{Q}{}^{({K_{s}})})$ satisfies

Theorem 3.2 states that for appropriately chosen threshold levels and all thresholding function satisfying (4), after enough iterations, Algorithm 1 yields principal subspace estimators whose errors are, with high probability, uniformly bounded over $\mathcal{F}_{n}$ by a sequence of asymptotically vanishing constants as $n\to\infty$. In addition, the probability that the estimation error is not well controlled vanishes polynomially fast. Therefore, the subspace estimators are uniformly consistent over $\mathcal{F}_{n}$.

The interpretation of the two terms in the error bound (16) is similar to those in Theorem 3.1. Having introduced those quantities in Section 3.3, we could elaborate a little more on the first, that is, the “nonparametric” term. By Theorem 3.3 below, when estimating $\mathcal{P}_{m}$, Algorithm 1 focuses only on the coordinates in $H$, whose cardinality is $\operatorname{card}(H)=O(M_{n})$. Though $H$ does not appear explicitly in the rates, the rates depend crucially on its cardinality which is further upper bounded by $M_{n}$. Since $\tau_{nm}^{2}$ can be interpreted as the average error per coordinate, the total estimation error accumulated over all coordinates in $H$ is thus of order $O(M_{n}\tau_{nm}^{2})$. Moreover, as we will show later, the squared bias induced by focusing only on $H$ is also of order $O(M_{n}\tau_{nm}^{2})$. Thus, this term indeed comes from the bias-variance tradeoff of the nonparametric estimation procedure. The meaning of the second, that is, the “parametric,” term is the same as in Theorem 3.1. Finally, we note that both terms vanish as $n\to\infty$ under Conditions GR, SP and $\operatorname{AD}(m,\kappa)$.

The threshold levels $\alpha_{n}$ and $\gamma_{nj}$ in (7) and (8) as well as $K_{s}$ in (9) do not depend on unknown parameters. So the estimation procedure achieves the rates adaptively over a wide range of high-dimensional sparse settings.

In addition, (15) implies that Algorithm 1 only needs a relatively small number of iterations to yield the desired estimator. In particular, when the largest spike $\lambda_{1}^{2}$ is bounded away from zero, (15) shows that it suffices to have $K_{s}\asymp\log{n}$ iterations. We remark that it is not critical to run precisely $K_{s}$ iterations. The result holds when we stop anywhere between $K$ and $2K$.

Theorem 3.2 could also be extended to an upper bound for the risk. Note that $p_{n}^{-2}=o(\tau_{nm}^{2}\vee\varepsilon_{nm}^{2})$, and that the loss function (3) is always bounded above by $1$. The following result is a direct consequence of Theorem 3.2.

Correct exclusion property

We now switch to the model selection property of Algorithm 1. By the discussion in Section 2, an important motivation for the iterative thresholding procedure is to trade bias for variance by keeping low signal coordinates out of the orthogonal iterations. More specifically, it is desirable to restrict our effort to estimating those coordinates in $H$ and simply estimating those coordinates in $L$ with zeros.

By construction, Algorithm 2 yields an initial matrix with a lot of zeros, but Algorithm 1 is at liberty to introduce new nonzero coordinates. The following result shows that with high probability all the nonzero coordinates introduced are in the set $H$.

Under the setup of Theorem 3.2, uniformly over $\mathcal{F}_{n}$, with probability at least $1-C_{0}p_{n}^{-2}$, for all $k=0,\ldots,{K_{s}}$, the orthonormal matrix $\widehat{Q}{}^{(k)}$ has zeros in all its rows indexed by $L$, that is, $\widehat{Q}{}^{(k)}_{L\cdot}=0$.

We call the property in Theorem 3.3 “correct exclusion,” because it ensures that all the low signal coordinates in $L$ are correctly excluded from iterations. In addition, Theorem 3.3 shows that the principal subspace estimator is indeed spanned by a set of sparse loading vectors, where all loadings in $L$ are exactly zero.

Note that the initial matrix $\widehat{Q}{}^{(0)}$ has all its nonzero coordinates in $B$, which, with high probability, only selects “big” coefficients in the leading eigenvectors, whose magnitudes are no less than $O([\log{p_{n}}/(n\lambda_{m}^{4})]^{1/4})$. On the other hand, the set $H$ includes all coordinates with magnitude no less than $O([\log p_{n}/(nh(\lambda_{m}^{2}))]^{1/2})$. Thus, the minimum signal strength for $H$ is of smaller order than that for $B$. So, with high probability, $B$ is a subset of $H$ consisting only of its coordinates with “big” signals. Thus, though $\widehat{Q}{}^{(0)}$ excludes all the coordinates in $L$, it only includes “big” coordinates in $H$ and fails to pick those medium sized ones which are crucial for obtaining the convergence rate (16). Algorithm 1 helps to include more coordinates in $H$ along iterations and hence achieves (16).

Rates of convergence for individual eigenvector estimation

The primary focus of this paper is on estimating principal subspaces. However, when an individual eigenvector, say $q_{j}$, is identifiable, it is also of interest to see whether Algorithm 1 can estimate it well. The following result shows that for $K_{s}$ in (9), the $j$th column of $\widehat{Q}{}^{({K_{s}})}$ estimates $q_{j}$ well, provided that the $j$th eigenvalue is well separated from the rest of the spectrum.

Under the setup of Theorem 3.2, suppose for some $j\leq m$, both conditions $\operatorname{AD}(j-1,\kappa^{\prime})$ and $\operatorname{AD}(j,\kappa^{\prime})$ are satisfied for some constant $\kappa^{\prime}<\lim_{n\to\infty}\lambda_{1}^{2}/(\lambda_{m}^{2}-\lambda_{m+1}^{2})$. Then uniformly over $\mathcal{F}_{n}$, with probability at least $1-C_{0}p_{n}^{-2}$, $\widehat{q}_{j}^{({K_{s}})}$, the $j$th column of $\widehat{Q}{}^{({K_{s}})}$, satisfies

5 Choice of m𝑚m

The main results in this section are stated with the assumption that the subspace dimension $m$ is given. In what follows, we discuss how to choose $m$ and also how to estimate $\bar{m}$ based on data.

be an estimator for $\bar{m}$, where for any positive integer $k$,

Then in Algorithm 1, for a large constant $\bar{\kappa}$, we define

For $\widehat{\bar{m}}$ in (17), we have the following results.

Suppose Conditions GR and SP hold. Let $\widehat{\bar{m}}$ be defined in (17) with $B$ obtained by Algorithm 2 with $\alpha_{n}$ specified by (7) for some $\alpha>2\sqrt{3}$. Then, uniformly over $\mathcal{F}_{n}$, with probability at least $1-C_{0}p_{n}^{-2}$: {longlist}[(3)]

for any $m$ such that the condition $\operatorname{AD}(m,\kappa)$ is satisfied with some constant $\kappa$, $m\leq\widehat{\bar{m}}$ when $n$ is sufficiently large;

if the condition $\operatorname{AD}(\bar{m},\kappa)$ is satisfied with some constant $\kappa$, $\widehat{\bar{m}}=\bar{m}$ when $n$ is sufficiently large.

By claim (1), any $m\leq\widehat{\bar{m}}$ satisfies $m\leq\bar{m}$ with high probability. In addition, claim (2) shows that, for sufficiently large $n$, any $m$ such that $\operatorname{AD}(m,\kappa)$ holds is no greater than $\widehat{\bar{m}}$. Thus, when restricting to those $m\leq\widehat{\bar{m}}$, we do not miss any $m$ such that Theorems 3.2 and 3.3 hold for estimating $\mathcal{P}_{m}$. These two claims jointly ensure that we do not need to consider any target dimension beyond $\widehat{\bar{m}}$. Finally, claim (3) shows that we recover the exact number of spikes with high probability for large samples when $\operatorname{AD}(\bar{m},\kappa)$ is satisfied, that is, when $\lambda_{1}^{2}\asymp\lambda_{\bar{m}}^{2}$. Note that this assumption was made in pajo07 .

Computational complexity

We now study the computational complexity of Algorithm 1. Throughout, we assume the same setup as in Section 3, and restrict the calculation to the high probability event on which the conclusions of Theorems 3.2 and 3.3 hold. For any matrix $A$, we use $\operatorname{supp}\{A\}$ to denote the index set of the nonzero rows of $A$.

Consider a single iteration, say, the $k$th. In the multiplication step, the $(\nu,j)$th element of $T^{(k)}$, $t_{\nu j}^{(k)}$, comes from the inner product of the $\nu$th row of $S$ and the $j$th column of $\widehat{Q}{}^{(k-1)}$. Though both are $p$-vectors, Theorem 3.3 asserts that for any column of $\widehat{Q}{}^{(k-1)}$, at most $\operatorname{card}(H)$ of its entries are nonzero. So if we know $\operatorname{supp}\{\widehat{Q}{}^{(k-1)}\}$, then $t_{\nu j}^{(k)}$ can be calculated in $O(\operatorname{card}(H))$ flops, and $T^{(k)}$ in $O(mp\operatorname{card}(H))$ flops. Since $\operatorname{supp}\{\widehat{Q}{}^{(k-1)}\}$ can be obtained in $O(mp)$ flops, the multiplication step can be completed in $O(mp\operatorname{card}(H))$ flops. Next, the thresholding step performs elementwise operation on $T^{(k)}$, and hence can be completed in $O(mp)$ flops. Turn to the QR step. First, we can obtain $\operatorname{supp}\{\widehat{T}{}^{(k)}\}$ in $O(mp)$ flops. Then QR factorization can be performed on the reduced matrix which only includes the rows in $\operatorname{supp}\{\widehat{T}{}^{(k)}\}$. Since Theorem 3.3 implies $\operatorname{supp}\{\widehat{T}{}^{(k)}\}=\operatorname{supp}\{\widehat{Q}{}^{(k)}\}\subset H$, the complexity of this step is $O(m^{2}\operatorname{card}(H))$. Since $m=O(p)$, the complexity of the multiplication step dominates, and so the complexity of each iteration is $O(mp\operatorname{card}(H))$. Theorem 3.2 shows that $K_{s}$ iteration is enough. Therefore, the overall complexity of Algorithm 1 is $O({K_{s}}mp\operatorname{card}(H))$.

When the true eigenvectors are sparse, $\operatorname{card}(H)$ is of manageable size. In many realistic situations, $\lambda_{1}^{2}$ is bounded away from and so $K_{s}\asymp\log{n}$. For these cases, Algorithm 1 is scalable to very high dimensions.

We conclude the section with a brief discussion on parallel implementation of Algorithm 1. In the $k$th iteration, both matrix multiplication and elementwise thresholding can be computed in parallel. For QR factorization, one needs only to communicate the rows of $\widehat{T}{}^{(k)}$ with nonzero elements, the number of which is no greater than $\operatorname{card}(H)$. Thus, the overhead from communication is $O(m\operatorname{card}(H))$ for each iteration, and $O({K_{s}}m\operatorname{card}(H))$ in total. When the leading eigenvectors are sparse, $\operatorname{card}(H)$ is manageable, and parallel computing of Algorithm 1 is feasible.

Numerical experiments

We first consider the case where each $x_{i}$ is generated by (6) with ${\bar{m}}=1$. Motivated by functional data with localized features, four test vectors $q_{1}$ are considered, where $q_{1}=(f(1/p),\ldots,f(p/p))^{\prime}$, with $f$ one of the four functions in Figure 1. For each test vector, the dimension $p=2048$, the sample size $n=1024$ and $\lambda_{1}^{2}$ ranges in $\{100,25,10,5,2\}$.

Before applying any sparse PCA method, we transform the observed data vectors into the wavelet domain using the Symmlet 8 basis mall09 , and scale all the observations by $\widehat{\sigma}$ with $\widehat{\sigma}^{2}$ given in (5). The multi-resolution plots of wavelet coefficients of the test vectors are shown in Figure 2. In the wavelet domain, the four vectors exhibits different levels of sparsity, with step the least sparse, and sing the most.

Table 1 compares the average loss of subspace estimation over $100$ runs for each spike value and each test vector by Algorithm 1 (ITSPCA) with several existing methods: augmented sparse PCA (AUGSPCA) pajo07 , correlation augmented sparse PCA (CORSPCA) nadl09 and diagonal thresholding sparse PCA (DTSPCA) given in Algorithm 2. For ITSPCA, we computed $\widehat{Q}{}^{(0)}$ by Algorithm 2. $\alpha_{n}$ and $\gamma_{n1}$ are specified by (7) and (8) with $\alpha=3$ and $\gamma=1.5$. These values are smaller than those in theoretical results, but lead to better numerical performance. We stop iterating once $L(\operatorname{ran}(\widehat{Q}{}^{(k)}),\operatorname{ran}(\widehat{Q}{}^{(k+1)}))\leq n^{-2}$. Parameters in competing algorithms are all set to the values recommended by their authors.

From Table 1, ITSPCA and CORSPCA outperform the other two methods in all settings. Between the two, CORSPCA only wins by small margins when the spike values are large. Otherwise, ITSPCA wins, sometimes with large margins. For the same algorithm at the same spike value, the sparser the signal, the smaller the estimation error.

Table 1 also presents the average sizes of the sets of selected coordinates. While all methods yield sparse PC loadings, AUGSPCA and DTSPCA seem to select too few coordinates, and thus introduce too much bias. ITSPCA and CORSPCA apparently result in a better bias-variance tradeoff.

2 Multiple spike settings

Next, we simulated data vectors using model (6) with ${\bar{m}}=4$. The $q_{j}$ vectors are taken to be the four test vectors used in single spike settings, in the same order as in Figure 1, up to orthonormalization.The four test vectors are shifted such that the inner product of any pair is close to . So the vectors after orthonormalization are visually indistinguishable from those in Figure 1. We tried four different configurations of the spike values $(\lambda_{1}^{2},\ldots,\lambda_{4}^{2})$, as specified in the first column of Table 2. For each configuration of spike values, the dimension is $p=2048$, and the sample size is $n=1024$.

For each simulated dataset, we estimate $\mathcal{P}_{m}$ for $m=1,2,3$ and $4$. The last four columns of Table 2 present the losses in estimating subspaces, averaged over $100$ runs, using the same sparse PCA methods as in single spike settings. For ITSPCA, we set the thresholds $\{\gamma_{nj},j=1,\ldots,4\}$ as in (8) with $\gamma=1.5$. All other implementation details are the same. Again, we used recommended values for parameters in all other competing methods.

The simulation results reveal two interesting phenomena. First, when the spikes are relatively well separated (the first and the last blocks of Table 2), all methods yield decent estimators of $\mathcal{P}_{m}$ for all values of $m$, which implies that the individual eigenvectors are also estimated well. In this case, ITSPCA always outperforms the other three competing methods. Second, when the spikes are not so well separated (the middle two blocks, with $m=1,2$ or $3$), no method leads to decent subspace estimator. However, all methods give reasonable estimators for $\mathcal{P}_{4}$ because $\lambda_{4}^{2}$ in both cases are well above . This implies that, under such settings, we fail to recover individual eigenvectors, but we can still estimate $\mathcal{P}_{4}$ well. ITSPCA again gives the smallest average losses. In all configurations, the estimated number of spikes $\widehat{\bar{m}}$ in (17) and the data-based choice of $m$ in (20) with $\bar{\kappa}=15$ consistently picked $m=\widehat{\bar{m}}=4$ in all simulated datasets. Therefore, we are always led to estimating the “right” subspace $\mathcal{P}_{4}$, and ITSPCA performs favorably over the competing methods.

In summary, simulations under multiple spike settings not only demonstrate the competitiveness of Algorithm 1, but also suggest: {longlist}[(2)]

The quality of principal subspace estimation depends on the gap between successive eigenvalues, in addition to the sparsity of eigenvectors;

Focusing on individual eigenvectors can be misleading for the purpose of finding low-dimensional projections.

Proof

This section is devoted to the proofs of Theorems 3.2 and 3.3. We state the main ideas in Section 6.1 and divide the proof into three major steps, which are then completed in sequel in Sections 6.2–6.4. Others results in Section 3.4 are proved in the supplementary material supp .

Major steps of the proof

In the $k$th iteration of the oracle Algorithm 1, denote the matrices obtained after multiplication and thresholding by

A joint proof of Theorems 3.2 and 3.3 can then be completed by the following three major steps: {longlist}[(3)]

In what follows, we complete the three steps in Sections 6.2–6.4.

Consider the “bias” part first. Define the oracle covariance matrix

A proof is given in the supplementary material supp . Weyl’s theorem (stsu90 , Corollary 4.4.10) and Davis–Kahn’s $\sin\theta$ theorem daka70 are the key ingredients in the proof here, and also in the proofs of Lemmas 6.2 and 6.3. Here, claim (1) only requires Conditions GR and SP, but not the condition $\operatorname{AD}(m,\kappa)$.

3 Properties of the oracle sequence

Uniformly over $\mathcal{F}_{n}$, with probability at least $1-C_{0}p_{n}^{-2}$: {longlist}[(4)]

for sufficiently large $n$, $K_{s}\in[K,2K]$.

A proof is given in the supplementary material supp . Here, claims (1) and (2) do not require the condition $\operatorname{AD}(m,\kappa)$. In claim (3), the bound $(1-\rho)^{2}/5$ is much larger than that in (16). For instance, if $\lambda_{m}^{2}+1\asymp\lambda_{m}^{2}-\lambda_{m+1}^{2}$, Lemmas 6.1 and 6.2 imply that $(1-\rho^{2})/5\asymp 1$ with high probability.

Claims (1) and (2) here, together with claims (1) of Lemmas 6.1 and 6.2, lead to the following result on consistent estimation of $\lambda_{1}^{2}/(\lambda_{j}^{2}-\lambda_{j+1}^{2})$ and $\lambda_{j}^{2}$, the proof of which is given in the supplementary material supp .

Evolution of the oracle sequence

The following proposition describes the evolution of ${\theta}^{(k)}$ over iterations.

Let $n$ be sufficiently large. On the event such that the conclusions of Lemmas 6.1–6.3 hold, uniformly over $\mathcal{F}_{n}$, for all $k\geq 1$: {longlist}[(2)]

then so is $\sin^{2}{{\theta}^{(k)}}$. Otherwise,

Together with Lemma 6.3, Proposition 6.1 also justifies the previous claim that elements of the oracle sequence are orthonormal with high probability.

Convergence

Finally, we study how fast the oracle sequence converges to a stable subspace estimator, and how good this estimator is.

For sufficiently large $n$, on the event such that the conclusions of Lemmas 6.1–6.3 hold, uniformly over $\mathcal{F}_{n}$, it takes at most $K$ steps for the oracle sequence to converge. In addition, there exist constants $C_{1}=C_{1}(\gamma,r,m,\kappa)$ and $C_{2}$, such that for all $k\geq K$,

A proof is given in the supplementary material supp , and this completes step 2.

4 Proof of main results

We now prove the properties of the actual estimating sequence. The proof relies on the following lemma, which shows the actual and the oracle sequences are identical up to $2K$ iterations.

A proof is given in the supplementary material supp , and this completes step 3.

We now prove Theorems 3.2 and 3.3 by showing that the actual sequence inherits the desired properties from the oracle sequence. Since Theorem 3.1 is a special case of Theorem 3.2, we do not give a separate proof.

Proof of Theorem 3.2 Note that the event on which the conclusions of Lemmas 6.1–6.4 hold has probability at least $1-C_{0}p_{n}^{-2}$. On this event,

Here, the first equality comes from Lemma 6.4. The first two inequalities result from the triangle inequality and Jensen’s inequality, respectively. Finally, the last inequality is obtained by noting that $K_{s}\in[K,2K]$ and by replacing all the error terms by their corresponding bounds in Lemmas 6.1, 6.2 and Proposition 6.2.

Acknowledgment

The author would like to thank Iain Johnstone for many helpful discussions.

Supplement to “Sparse principal component analysis and iterative thresholding” \slink[doi]10.1214/13-AOS1097SUPP \sdatatype.pdf \sfilenameaos1097_supp.pdf \sdescriptionWe give in the supplement proofs to Corollaries 3.1 and 3.2, Proposition 3.1 and all the claims in Section 6.

References