Sparse Principal Components Analysis

Introduction

Suppose $\{x_{i},i=1,\ldots,n\}$ is a dataset of $n$ observations on $p$ variables. Standard principal components analysis (PCA) looks for vectors $\xi$ that maximize

If $\xi_{1},\ldots,\xi_{k}$ have already been found by this optimization, then the maximum defining $\xi_{k+1}$ is taken over vectors $\xi$ orthogonal to $\xi_{1},\ldots,\xi_{k}$.

Our interest lies in situations in which each $x_{i}$ is a realization of a possibly high dimensional signal, so that $p$ is comparable in magnitude to $n$, or may even be larger. In addition, we have in mind settings in which the signals $x_{i}$ contain localized features, so that the principal modes of variation sought by PCA may well be localized also.

Consider, for example, the sample of an electrocardiogram (ECG) in Figure 1 showing some 13 consecutive heart beat cycles as recorded by one of the standard ECG electrodes. Individual beats are notable for features such as the sharp spike (“QRS complex”) and the subsequent lower peak (“T wave”), shown schematically in the second panel. The presence of these local features, of differing spatial scales, suggests the use of wavelet bases for efficient representation. Traditional ECG analysis focuses on averages of a series of beats. If one were to look instead at beat to beat variation, one might expect these local features to play a significant role in the principal component eigenvectors.

Returning to the general situation, the main contentions of this paper are:

(a) that when $p$ is comparable to $n$, some reduction in dimensionality is desirable before applying any PCA-type search for principal modes, and

(b) the reduction in dimensionality is best achieved by working in a basis in which the signals have a sparse representation.

We will support these assertions with arguments based on statistical performance and computational cost.

We begin, however, with an illustration of our results on a simple constructed example. Consider a single component (or single factor) model, in which, when viewed as $p-$dimensional column vectors

Panel (a) of Figure 2 shows an example of $\rho$ with $p=2048$ and the vector $\rho_{l}=f(l/n)$ where $f(t)$ is a mixture of Beta densities on $$, scaled so that$\|\rho\|=(\sum_{1}^{p}\rho_{l}^{2})^{1/2}=10.$Panels (b) and (c) show two sample paths from model (2): the random effect$v_{i}\rho$is hard to discern individual cases. Panel (d) shows the result of standard PCA applied to$n=1024$observations from (2) with$\sigma=1$. The effect of the noise remains clearly visible in the estimated principal eigenvector.

For functional data of this type, a regularized approach to PCA has been proposed by Rice & Silverman (1991) and Silverman (1996), see also Ramsay & Silverman (1997) and references therein. While smoothing can be incorporated in various ways, we illustrate the method discussed also in Ramsay & Silverman (1997, Ch. 7), which replaces (1) with

where $D^{2}\xi$ is the $(p-2)\times 1$ vector of second differences of $\xi$ and $\lambda\in(0,\infty)$ is the regularization parameter.

Panel (e) shows the estimated first principal component vector found by maximizing (3) with $\lambda=10^{-12}$ and $\lambda=10^{-6}$ respectively. Neither is really satisfactory as an estimate: the first recovers the original peak heights, but fails fully to suppress the remaining baseline noise, while the second grossly oversmooths the peaks in an effort to remove all trace of noise. Further investigation with other choices of $\lambda$ confirms the impression already conveyed here: no single choice of $\lambda$ succeeds both in preserving peak heights and in removing baseline noise.

Panel (f) shows the result of the adaptive sparse PCA algorithm to be introduced below: evidently both goals are accomplished quite satisfactorily in this example.

The need to select subsets: (in)consistency of classical PCA

A basic element of our sparse PCA proposal is initial selection of a relatively small subset of the initial $p$ variables before any PCA is attempted. In this section, we formulate some (in)consistency results that motivate this initial step.

Consider first the single component model (2). The presence of noise means that the sample covariance matrix $S=n^{-1}\sum_{i=1}^{n}x_{i}x_{i}^{T}$ will typically have $\min(n,p)$ non-zero eigenvalues. Let $\hat{\rho}$ be the unit eigenvector associated with the largest sample eigenvalue—with probability one it is uniquely determined up to sign.

One natural measure of the closeness of $\hat{\rho}$ to $\rho$ uses the angle $\angle(\hat{\rho},\rho)$ between the two vectors. We decree that the signs of $\hat{\rho}$ and $\rho$ be taken so that $\angle(\hat{\rho},\rho)$ lies in $[0,\pi/2]$. It will be convenient to phrase the results in terms of an equivalent distance measure

For asymptotic results, we will assume that there is a sequence of models (2) indexed by $n$. Thus, we allow $p(n)$ and $\rho(n)$ to depend by $n$, though the dependence will usually not be shown explicitly. [Of course $\sigma$ might also be allowed to vary with $n$, but for simplicity it is assumed fixed.]

Our first interest is whether the estimate $\hat{\rho}$ is consistent as $n\rightarrow\infty$. This turns out to depend crucially on the limiting value

One setting in which this last assumption may be reasonable is when $p(n)$ grows by adding finer scale wavelet coefficients of a fixed function as $n$ increases.

Then with probability one as $n\rightarrow\infty$,

so long as the right side is at most one.

For the proof, see Appendix A.2. The bound $\zeta(\tau;c)$ is decreasing in the “signal-to-noise” ratio $\tau=\varrho/\sigma$ and increasing in the dimension-to-sample size ratio $c=\lim p/n$. It approaches as $c\rightarrow 0$, and in particular it follows that $\hat{\rho}$ is consistent if $p/n\rightarrow 0$.

The proof is based on an almost sure bound for eigenvectors of perturbed symmetric matrices. It appears to give the correct order of convergence: in the case $p/n\rightarrow 0$, we have

with $c(\tau)=4\tau^{-1}+8\tau^{-2}$, and examination of the proof shows that in fact

which is consistent with the $n^{-1/2}$ convergence rate that is typical when $p$ is fixed.

However if $c>0$, the upper bound (7) is strictly positive. And it turns out that $\hat{\rho}$ must be an inconsistent estimate in this setting:

Assume model (2), (5) and (6). If $p/n\rightarrow c>0$, then $\hat{\rho}$ is inconsistent:

In short, $\hat{\rho}$ is a consistent estimate of $\rho$ if and only if $p=o(n)$. The noise does not average out if there are too many dimensions $p$ relative to sample size $n$. A heuristic explanation for this phenomenon is given just before the proof in Appendix A.3.

The inconsistency criterion extends to a considerably more general multi-component model. Assume that we have $n$ curves $x_{i}$, observed at $p$ time points. Viewed as $p$ dimensional column vectors, this model assumes that

Here $\mu$ is the mean function, which is assumed known, and hence is taken to be zero. We make the following assumptions:

(a) The $\rho^{j},j=1,...,m\leq p$ are unknown, mutually orthogonal principal components, with norms $\rho_{j}(n)=\|\rho^{j}\|$

(b) The multipliers $v_{i}^{j}\sim N(0,1)$ are all independent over $j=1,\ldots,m$ and $i=1,\ldots,m$.

(c) The noise vectors $z_{i}\sim N_{p}(0,I)$ are independent among themselves and also of the random effects $\{v_{i}^{j}\}$.

We continue to focus on the estimation of the principal eigenvector $\rho^{1}$, and establish a more general version of the two preceding theorems.

Assume model (8) together with conditions (a)-(d). If $p/n\rightarrow c$, then

so long as the right side is at most, say, 4/5.

Thus, it continues to be true in the multicomponent model that $\hat{\rho}^{1}$ is consistent if and only if $p=o(n)$.

The sparse PCA algorithm

The inconsistency results of Theorems 2 and 3 emphasize the importance of reducing the number of variables before embarking on PCA, and motivate the sparse PCA algorithm to be described in general terms here. Note that the algorithm per se does not require the specification of a particular model, such as (8).

[The wavelet basis is used in this paper, for reasons discussed in the next subsection.]

2. Subset. Calculate the sample variances $\hat{\sigma}_{\nu}^{2}=\widehat{Var}(x_{i\nu})$. Let $\hat{I}$ denote the set of indices $\nu$ corresponding to the largest $k$ variances.

[$k$ may be specified in advance, or chosen based on the data, see Section 3.2 below].

3. Reduced PCA. Apply standard PCA to the reduced data set $\{x_{i\nu},\nu\in\hat{I},i=1,\ldots,n\}$ on the selected $k-$dimensional subset, obtaining eigenvectors $\hat{\rho}^{j}=(\hat{\rho}^{j}_{\nu}),j=1,\ldots,k$.

4. Thresholding. Filter out noise in the estimated eigenvectors by hard thresholding

[Hard thresholding is given, as usual, by $\eta_{H}(x,\delta)=xI\{|x|\geq\delta\}$. An alternative is soft thresholding $\eta_{S}(x,\delta)=\text{sgn}(x)(|x|-\delta)_{+}$, but hard thresholding has been used here because it preserves the magnitude of retained signals.

The threshold $\delta$ can be chosen, for example, by trial and error, or as $\delta=\hat{\tau}_{j}\sqrt{2\log k}$ for some estimate $\hat{\tau}_{j}.$ In this paper, estimate (13) is used. Another possibility is to set $\hat{\tau}_{j}=MAD\{\hat{\rho}^{j}_{\nu},\nu=1,\ldots,k\}/0.6745$. ]

5. Reconstruction. Return to the original signal domain, setting

In the rest of this section, we amplify on and illustrate various aspects of this algorithm. Given appropriate eigenvalue and eigenvector routines, it is not difficult to code. For example, MATLAB files that produce most figures in this paper will soon be available at www-stat.stanford.edu/~imj/ – to exploit wavelet bases, they make use of the open-source library WaveLab available at www-stat.stanford.edu/~wavelab/.

Suppose that in the basis $\{e_{\nu}(t)\}$ a population principal component $\rho(t)$ has coefficients $\{\rho_{\nu}\}$:

It is desirable, both from the point of view of economy of representation, as well as computational complexity, for the expansion in basis $\{e_{\nu}\}$ to be sparse, i.e., most coefficients $\rho_{\nu}$ are small or zero.

Wavelet bases typically provide sparse representations of one-dimensional functions that are smooth or have isolated singularities or transient features, such as in our ECG example. Here is one such result. Expand $\rho$ in a nice wavelet basis $\{\psi_{jk}(t)\}$ to obtain $\rho=\sum_{jk}\rho_{jk}\psi_{jk}(t)$ and then order coefficients by absolute magnitude, so that $(\rho_{\nu})$ is a re-ordering of the $|\rho_{jk}|$ in decreasing order. Then smoothness (as measured by membership in some Besov space $B^{\alpha}_{p,q}$) implies sparsity in the sense that

[for details, see Donoho (1993) and Johnstone (2002): in particular it is assumed that $\alpha>(1/p-1/2)_{+}$ and that the wavelet $\psi$ is sufficiently smooth.]

In this paper, we will assume that the basis $\{e_{\nu}\}$ is fixed in advance – and it will generally be taken to be a wavelet basis. Extension of our results to incorporate basis selection (e.g. from a library of orthonormal bases such as wavelet packets) is a natural topic for further research.

2 Adaptive choice of k𝑘k

Here are two possibilities for adaptive choice of $\hat{k}=|\hat{I}|$ from the data:

(a) choose co-ordinates with variance exceeding the estimated noise level by a specified fraction $\alpha_{n}$:

This choice is considered further in Section 3.5.

(b) As motivation, recall that we hope that the selected set of variables $\hat{I}$ is both small in cardinality and also captures most of the variance of the population principal components, in the sense that the ratio

is close to one for the leading population principal components in $\{\rho^{1},\ldots,\rho^{m}\}$. Now let $\chi_{(n),\alpha}^{2}$ denote the upper $\alpha-$percentile of the $\chi_{(n)}^{2}$ distribution – if all co-ordinates were pure noise, one might expect $\hat{\sigma}_{(\nu)}^{2}$ to be close to $n^{-1}\hat{\sigma}^{2}\chi_{(n),\nu/n}^{2}$. Define the excess over these percentiles by

where $\hat{k}$ is the smallest index $k$ for which the inequality holds. This second method has been used for the figures in this paper, typically with $w(n)=.995$.

Estimation of $\sigma$. If the population principal components $\rho^{j}$ have a sparse representation in basis $\{e_{\nu}\}$, then we may expect that in most co-ordinates $\nu$, $\{x_{i\nu}\}$ will consist largely of noise. This suggests a simple estimate of the noise level on the assumption that the noise level is the same in all co-ordinates, namely

3 Computational complexity

It is straightforward to estimate the cost of sparse PCA by examining its main steps:

This depends on the choice of basis. In the wavelet case no more than $O(np\log p)$ operations are needed.

Sort the sample variances and select $\hat{I}$: $O(p\log p)$.

Eigendecomposition for a $k\times k$ matrix: $O(k^{3})$.

Estimate $\hat{\sigma}^{2}$ and $\widehat{\|\rho\|^{2}}$: $O(p)$.

Reconstruct eigenvectors in the original sample space: $O(k^{2}p)$.

Both standard and smoothed PCA need at least $O((p\wedge n)^{3})$ operations. Therefore, if we can find a sparse basis such that $k/p\rightarrow 0$, then under the assumption that $p/n\rightarrow c$ as $n\rightarrow\infty$,the total cost of sparse PCA is $o(p^{3})$. We will see in examples to follow that the savings can be substantial.

4 Simulated examples

The two examples in this section are both motivated by functional data with localized features.

The first is a three-peak principal component depicted in Figure 2, and already discussed in Section 1. The second example, Figure 3, has an underlying first principal component composed of step functions. For both examples, the dimension of data vectors is $p=2048$, the number of observations $n=1024$, and the noise level $\sigma=1$. However, the amplitudes of $\rho$ differ, with $\|\rho\|=10$ for the “3-peak” function and $\|\rho\|\approx 25$ for the “step” function.

Panels (d) and (b) in the two figures respectively show the sample principal components obtained by using standard PCA. While standard PCA does capture the peaks and steps, it retains significant noise in the flat regions of the function. Corresponding panels (e) and (c) show results from smooth PCA with the indicated values of the smoothing parameter. Just as for the three peak curve discussed earlier, in the case of the step function, none of the three estimates simultaneously captures both jumps and flat regions well.

Panels (f) and (d) present the principal components obtained by sparse PCA. Using method (b) of the previous section with $w=99.5\%$, the Subset step selects $k=372$ and 438 for the “3-peak” curve and “step” function, respectively. The sample principal component in Figure 2(d) is clearly superior to the other sample p.c.s in Figure 2. Although the principal component function in the step case appears to be only slightly better than the solid blue smooth PCA estimate, we will see later that its squared error is reduced by more than 90%.

Table 1 compares the accuracy of the three PCA algorithms, using average squared error (ASE) defined as

The average ASE over 50 iterations is shown. The running time is the CPU time for a single iteration used by Matlab on a MIPS R10000 195.0MHz server.

Figure 4 presents box plots of ASE for the 50 iterations. Sparse PCA gives the best result for the “step” curve. For the “3-peak” function, in only a few iterations does sparse PCA generate larger error than smoothed PCA with a small $\lambda=10^{-12}$. On the average, ASE using sparse PCA is superior to the other methods by a large margin. Overall Table 1 and Figure 5 show that sparse PCA leads to the most accurate principal component while using much less CPU time than other PCA algorithms.

Anderson (1963) obtained the asymptotic distribution of $\sqrt{n}(\rho-\hat{\rho})$ for fixed $p$; in particular

as $n\rightarrow\infty.$ For us, $p$ increases with $n$, but we will nevertheless use (12) as an heuristic basis for estimating the variance $\hat{\tau}$ needed for thresholding. Since the effect of thresholding is to remove noise in small coefficients, setting $\rho_{\nu}$ to 0 in $(\ref{eq:vardiff})$ suggests

Neither $\|\rho\|^{2}$ and $\sigma^{2}$ in $(\ref{eq:taui})$ are known, but they can be estimated by using the information contained in the sample covariance matrix $S$, much as in the discussion of Section 3.2. Indeed $S_{\nu}^{2}$, the $\nu$-th diagonal element of $S$, follows a scaled $\chi^{2}$ distribution, with expectation $\rho_{\nu}^{2}+\sigma^{2}.$ If $\rho_{\nu}$ is a sparse representation of $\rho$, then most coefficients will be small, suggesting the estimate (11) for $\sigma^{2}$. In the single component model,

Figure 5 shows the histograms for these estimates of $\|\rho\|$ and $\sigma$ based on 100 iterations for the “3-peak” curve and for the “step” function.

5 Correct Selection Properties

A basic issue raised by the sparse PCA algorithm is whether the selected subset $\hat{I}$ in fact correctly contains the largest population variances, and only those. We formulate a result, based on large deviations of $\chi^{2}$ variables, that provides some reassurance.

For this section, assume that the diagonal elements of the sample covariance matrix $S=n^{-1}\sum_{1}^{n}x_{i}x_{i}^{T}$ have marginal $\chi^{2}$ distributions, i.e.,

We will not require any assumptions on the joint distribution of $\{\hat{\sigma}_{\nu}^{2}\}$.

Denote the ordered population coordinate variances by $\sigma_{(1)}^{2}\geq\sigma_{(2)}^{2}\geq\ldots$ and the ordered sample coordinate variances by $\hat{\sigma}_{(1)}^{2}\geq\hat{\sigma}_{(2)}^{2}\geq\ldots$. A desirable property is that $\hat{I}$ should, for suitable $\alpha_{n}$ small,

We will show that this in fact occurs if $\alpha_{n}=\gamma\sqrt{n^{-1}\log n}$, for appropriate $\gamma>0.$

We say that a false exclusion (FE) occurs if any variable in $I_{in}$ is missed:

while a false inclusion (FI) happens if any variable in $I_{out}$ is spuriously selected:

Under assumptions (15), the chance of an inclusion error of either type in $\hat{I}_{k}$ having magnitude $\alpha_{n}=\gamma n^{-1/2}(\log n)^{1/2}$ is polynomially small:

with $b(\gamma)=[\gamma\sqrt{3}/(4+2\sqrt{3})]^{2}.$

For example, if $\gamma=9$, then $b(\gamma)\doteq 4.36.$ As a numerical illustration based on (54) below, if the subset size $k=50$, while $p=n=1000$, then the chance of an inclusion error corresponding to a 25% difference in SDs (i.e. $\sqrt{1+\alpha_{n}}=1.25$) is below 5%. That reasonably large sample sizes are needed is a sad fact inherent to variance estimation—as one of Tukey’s ‘anti-hubrisines’ puts it, “it takes 300 observations to estimate a variance to one significant digit of accuracy”.

6 Consistency

The sparse PCA algorithm is motivated by the idea that if the p.c.’s have a sparse representation in basis $\{e_{\nu}\}$, then selection of an appropriate subset of variables should overcome the inconsistency problem described by Theorem 2.

To show that such a hope is justified, we establish a consistency result for sparse PCA. For simplicity, we consider the single component model (2), and assume that $\sigma^{2}$ is known—though this latter assumption could be removed by estimating $\sigma^{2}$ using (11).

To select the subset of variables $\hat{I}$, we use a version of rule (a) from Section 3.2:

with $\gamma_{n}=\gamma(n^{-1}\log n)^{1/2}$ and $\gamma$ a sufficiently large positive constant—for example $\gamma>\sqrt{12}$ would work for the proof.

We assme that the unknown principal components $\rho=\rho(n)$ satisfy a uniform sparsity condition: for some positive constants $q,C$,

Let $\hat{\rho}_{I}$ denote the principal eigenvector estimated by step (3) of the sparse PCA algorithm (thresholding is not considered here).

Assume that the single component model (2) holds, with $p/n\rightarrow c>0$ and $\|\rho(n)\|\rightarrow\varrho>0$. For each $n$, assume that $\rho(n)$ satisfies the uniform sparsity condition (17).

Then the estimated principal eigenvector $\hat{\rho}_{I}$ obtained by subset selection rule (16) is consistent:

The proof is given in Appendix A.5: it is based on a correct selection property similar to Theorem 4: combined with a modification of the consistency argument for Theorem 3. In fact, the proof shows that consistency holds even under the weaker assumption $p=O(n^{a})$, for arbitrary $a>0$, so long as $\gamma=\gamma(a)$ is set sufficiently large.

7 ECG example

This section offers a brief illustration of sparse PCA as applied to some ECG data kindly provided by Jeffrey Froning and Victor Froelicher in the cardiology group at Palo Alto Veterans Affairs Hospital. Beat sequences – typically about 60 cycles in length – were obtained from some 15 normal patients: we have selected two for the preliminary illustrations here.

Data Preprocessing. Considerable preprocessing is routinely done on ECG signals before the beat averages are produced for physician use. Here we describe certain steps taken with our data, in collaboration with Jeff Froning, preparatory to the PCA analysis.

The most important feature of an ECG signal is the Q-R-S complex: the maximum occurs at the R-wave, as depicted in Figure 1(b). Therefore we define the length of one cycle as the gap between two adjacent maxima of R-waves.

1. Baseline wander is observed in many ECG data sets, c.f. Figure 6. One common remedy for this problem is to deduct a piecewise linear baseline from the signal, the linear segment (dashed line) between two beats being determined from two adjacent onset points.

The onset positions of R-waves are shown by asterisks. Their exact locations vary for different patients, and as Figure 6 shows, even for adjacent R-waves. The locations are determined manually in this example. To reduce the effect of noise, the values of onset points are calculated by an average of 5 points close to the onset position.

2. Since pulse rates vary even on short time scales, the duration of each heart beat cycle may vary as well. We use linear interpolation to equalize the duration of each cycle, and for convenience in using wavelet software, discretize to $512=2^{9}$ sample points in each cycle.

3. Finally, due to the importance of the R-wave, the horizontal positions of the maxima are the 150th position in each cycle.

4. Convert the ECG data vector into an $n\times p$ data matrix, where $n$ is the number of observed cycles and $p=512$. Each row of the matrix presents one heart beat cycle with the maxima of R-waves all aligned at the same position.

PCA analysis. Figure 7 (a) and (d) shows the mean curves for two ECG samples in blue. The number of observations $n$, i.e. number of heart beats recorded, are 66 and 61, respectively. The first sample principal components for these two sample sets are plotted in plots (c) and (f), with red curves from standard PCA and blue curves from sparse PCA. In both cases there are two sharp peaks in the vicinity of the QRS complex. The first peak occurs shortly before the 150th position, where all the maxima of R-waves are aligned, and the second peak, which has an opposite sign, shortly after.

The standard PCA curve in Figure 6.7.(b, red) is less noisy than that in panel (d, red), even allowing for the difference in vertical scales. Using $(\ref{eq:sigma2})$,

while the magnitudes of the two mean sample curves are very similar.

The sparse PCA curves (blue) are smoother than the standard PCA ones (red), especially in plot (d) where the signal to noise ratio is lower. On the other hand, the red and blue curves match quite well at the two main peaks. Sparse PCA has reduced noise in the sample principal component in the baseline while keeping the main features.

There is a notable difference between the two estimated p.c.’s. In the first case, the p.c. is concentrated around the R-wave maximum, and the effect is to accelerate or decelerate the rise (and fall) of this peak from baseline in a given cycle. This is more easily seen by comparing plots of $\bar{x}+2\hat{\rho}$ (green) with $\bar{x}-2\hat{\rho}$ (red), shown over a magnified part of the cycle in panel (b). In the second case, the bulk of the energy of the p.c. is concentrated in a level shift in the part of the cycle starting with the ST segment. This can be interpreted as beat to beat fluctuation in baseline – since each beat is anchored at at the onset point, there is less fluctuation on the left side of the peak. This is particularly evident in panel (e) – there is again a slight acceleration/deceleration in the rise to the R wave peak – less pronounced in the first case, and also less evident in the fall.

Obvious questions raised by this illustrative example include the nature of effects which may have been introduced by the preprocessing steps, notably the baseline removal anchored at onset points and the alignment of R-wave maxima. Clearly some conventions must be adopted to create rectangular data matrices for p.c. analysis, but detailed analysis of these issues must await future work.

Finally, sparse PCA uses less than 10% of the computing time than standard PCA.

Appendix A Appendix

Matrices. We first recall some pertinent matrix results. Define the $2-$norm of a rectangular matrix by

If $A$ is real and symmetric, then $\|A\|_{2}=\lambda_{max}(A).$ If $A_{p\times p}$ is partitioned

where $b$ is $(p-1)\times 1$, then by setting $x=(1\ 0^{T})^{T}$ in (18), one finds that

The matrix $B=\rho u^{T}+u\rho^{T}$ has at most two non-zero eigenvalues, given by

Indeed, the identity $\det(I+AC)=\det(I+CA)$ for compatible rectangular matrices $A$ and $C$ means that the non-zero eigenvalues of

are the same as those of the $2\times 2$ matrix

These remarks can be used to bound the angle between a vector $\eta$ and its image under a symmetric matrix $M$ in terms of the angle between $\eta$ and any principal eigenvector of $M$.

Let $\xi$ be a principal eigenvector of a non-zero symmetric matrix $M$. For any $\eta\neq 0$,

We may assume without loss of generality that $\|\xi\|=\|\eta\|=1$ and that $\xi^{T}\eta\geq 0$. Since $\xi$ is a principal eigenvector of a symmetric matrix, $\|M\xi\|=\|M\|$. From the sine rule (23),

where the final equality uses (22). Some calculus shows that $2\sin\alpha/2\leq\sin 2\alpha$ for $0\leq\alpha\leq\pi/4$ and hence

using (24) and the fact that $\xi$ is an eigenvector of $M$. ∎

Perturbation bounds. Suppose that a symmetric matrix $A_{p\times p}$ has unit eigenvector $q_{1}$. We wish to bound the effect of a symmetric perturbation $E_{p\times p}$ on $q_{1}$. The following result (Golub & Van Loan (1996, Thm 8.1.10), see also Stewart & Sun (1990)) constructs a unit eigenvector $\hat{q}_{1}$ of $A+E$ and bounds its distance from $q_{1}$ in terms of $\|E\|_{2}$. Here, the distance between unit eigenvectors $q_{1}$ and $\hat{q}_{1}$ is defined as at (4) and (21).

Let $Q_{p\times p}=[q_{1}\ Q_{2}]$ be an orthogonal matrix containing $q_{1}$ in the first column, and partition conformally

where $D_{22}$ and $E_{22}$ are both $(p-1)\times(p-1)$.

Suppose that $\lambda$ is separated from the rest of the spectrum of $A$; set

is a unit eigenvector of $A+E$. Moreover,

Let us remark that since $\|e\|_{2}\leq\|E\|_{2}$ by (19), we have $\|r\|_{2}\leq 1$ and

Suppose now that $q_{1}$ is the eigenvector of $A$ associated with the principal eigenvalue $\lambda_{1}(A)$. We verify that, under the preceding conditions, $\hat{q}_{1}$ is also the principal eigenvector of $A+E$: i.e. if $(A+E)\hat{q}_{1}=\lambda^{*}\hat{q}_{1}$, then in fact $\lambda^{*}=\lambda_{1}(A+E)$.

To show this, we verify that $\lambda^{*}>\lambda_{2}(A+E)$. Take inner products with $q_{1}$ in the eigenequation for $\hat{q}_{1}$:

Since $A$ is symmetric, $q_{1}^{T}A=\lambda_{1}(A)q_{1}^{T}$. Trivially, we have $q_{1}^{T}E\hat{q}_{1}\geq-\|E\|_{2}$. Combine these remarks with (26) to get

Now $\delta=\lambda_{1}(A)-\lambda_{2}(A)$ and since from the minimax characterization of eigenvalues (e.g. Golub & Van Loan (1996, p. 396) or Stewart & Sun (1990, p.218)), $\lambda_{2}(A+E)\leq\lambda_{2}(A)+\|E\|_{2}$, we have

Large Deviation Inequalities. If $\bar{X}=n^{-1}\sum_{1}^{n}X_{i}$ is the average of i.i.d. variates with moment generating function $\exp\{\Lambda(\lambda)\}=E\exp\{\lambda X_{1}\}$, then Cramer’s theorem (see e.g. Dembo & Zeitouni (1993, 2.2.2 and 2.2.12)) says that for $x>EX_{1}$,

where the conjugate function $\Lambda^{*}(x)=\sup_{\lambda}\{\lambda x-\Lambda(\lambda)\}$. The same bound holds for $P\{\bar{X}<x\}$ when $x<EX_{1}$.

When applied to the $\chi_{(n)}^{2}$ distribution, with $X_{1}=z_{1}^{2}$ and $z_{1}\sim N(0,1)$, the m.g.f. $\Lambda(\lambda)=-\textstyle{\frac{1}{2}}\log(1-2\lambda)$ and the conjugate function $\Lambda^{*}(x)=\textstyle{\frac{1}{2}}[x-1-\log x]$. The bounds

(the latter following, e.g., from (47) in Johnstone (2001)) yield

We will use also a slightly sharper bound

valid for $n\geq 16$ and $0\leq t\leq n^{1/6}$ (Johnstone, 2001).

When applied to sums of variables $X_{1}=z_{1}z_{2}$, with $z_{1}$ and $z_{2}$ independent $N(0,1)$ variates, the m.g.f. $\Lambda(\lambda)=-\textstyle{\frac{1}{2}}\log(1-\lambda^{2})$. With $\lambda_{*}(x)=[(1+4x^{2})^{1/2}-1]/(2x),$ the conjugate function satisfies

as $x\rightarrow 0$. Hence, for $n$ large,

Decomposition of sample covariance matrix. Now adopt the multicomponent model (8) along with its assumptions (a) - (c). The sample covariance matrix $S=n^{-1}\sum_{1}^{n}x_{i}x_{i}^{T}$ has expectation $ES=R+\sigma^{2}I_{p}$, where

Now decompose $S$ according to (8). Introduce $1\times n$ row vectors $v^{jT}=(v_{1}^{j}\cdots v_{n}^{j})$ and collect the noise vectors into a matrix $Z_{p\times n}=[z_{1}\cdots z_{n}]$. We then have

Some limit theorems. We turn to properties of the noise matrix $Z$ appearing in (35). The cross products matrix $ZZ^{T}$ has a standard $p$-dimensional Wishart $W_{p}(n,I)$ distribution with $n$ degrees of freedom and identity covariance matrix, see e.g. Muirhead (1982, p82). Thus the matrix $C=(c_{jk})$ in (34) is simply a scaled and recentered Wishart matrix. We state results below in terms of either $ZZ^{T}$ or $C$, depending on the subsequent application. Properties (b) and (c) especially play a key role in inconsistency when $c>0$.

We may clearly take $\sigma=1.$ An off-diagonal term in $n^{-1}ZZ^{T}=(c_{jk})$ has the distribution of an i.i.d. average $\bar{X}=n^{-1}\sum X_{i}$ where $X_{1}=z_{1}z_{2}$ is the product of two independent standard normal variates. Thus

Now apply the large deviation bound (32) to the right hand side. Since $p\sim cn$, the Borel-Cantelli lemma suffices to establish (36) for off-diagonal elements for any $b>2$.

A diagonal term $c_{jj}+1$ in $n^{-1}ZZ^{T}$ has the $n^{-1}\chi^{2}_{(n)}$ distribution. Setting $t=\sqrt{\textstyle{\frac{1}{2}}b\log n}$ in (31) yields

Since there are $p\sim cn$ diagonal terms, the conclusion (36) follows (again via Borel-Cantelli) so long as $b>8$. ∎

(b) Geman (1980) and Silverstein (1985) respectively established almost sure limits for the largest and smallest eigenvalues of a $W_{p}(n,I)$ matrix as $p/n\rightarrow c\in[0,\infty)$, from which follows:

[Although the results in the papers cited are for $c\in(0,\infty)$, the results are easily extended to $c=0$ by simple coupling arguments.]

(c) Suppose in addition that $v$ is a $1\times n$ vector with independent $N(0,1)$ entries, which are also independent of $Z$. Conditioned on $v$, the vector $Zv$ is distributed as $N_{p}(0,\|v\|^{2}I).$ Since $Z$ is independent of $v$, we conclude that

Now let $u_{p\times 1}=\sigma n^{-1}Zv$. From (39) we have

the distribution of the first component of $U_{p}$. It is well known that $U_{1}^{2}\sim\text{Beta}(1/2,(p-1)/2)$, so that $EU_{1}^{2}=p^{-1}$ and $\text{Var}U_{1}^{2}\leq 2p^{-2}$. From this it follows that

(d) Let $u^{j}=\sigma n^{-1}Zv^{j}$ be the vectors appearing in the definition of $B^{j}$ for $1\leq j\leq m$. We will show that a.s.

(the constant $c_{0}=2\sigma(1+\sqrt{c})$ would do).

From (38), it follows that w.p. 1, ultimately

The squared lengths $\|v^{j}\|^{2}$ follow independent $\chi_{(n)}^{2}$ laws. Since from (28) there exists $c_{1}$ for which $P\{\chi_{(n)}^{2}\geq 2n\}\leq e^{-c_{1}n}$ for $n\geq n_{0}$, it follows that

and so w.p. $1$ it is ultimately true that

Substituting (44) and (45) into (43), we recover (42). ∎

A.2 Upper Bounds: Proof of Theorems 1 and 3

Instead of working directly with the sample covariance matrix $S$, we consider $S^{*}=S-\sigma^{2}I_{p}.$ It is apparent that $S^{*}$ has the same eigenvectors as $S$. We decompose $S^{*}=R+E$, where $R$ is given by (33) and has spectrum

The perturbation matrix $E=A+B+C,$ where $A$ and $B$ refer to the sums in (34).

Assume that multicomponent model (8) holds, along with assumptions (a) - (d). For any $\epsilon>0$, if $p,n\rightarrow\infty,p/n\rightarrow c$, then almost surely

where the $A_{n},B_{n}$ and $C_{n}$ will be given below. We have shown explicitly the dependence on $\omega$ to emphasize that these quantities are random. Finally we show that the a.s. limit of $E_{n}(\omega)$ is the right side of (46).

The $v_{s}^{jk}$ are entries of a scaled and recentered $W_{m}(n,I)$ matrix, and so by (36), the maximum converges almost surely to . Since $\sum_{j}\|\rho^{j}\|\rightarrow\sum\varrho_{j}<\infty,$ it follows that the $A_{n}$-term converges to zero a.s.

$B$ term. Applying (20) to the definition (35) of $B^{j}$, we have

where $\tau_{j}=\rho^{jT}u^{j}/\|\rho^{j}\|\|u^{j}\|$ and $u^{j}=\sigma n^{-1}Zv^{j},$ and the convergence follows from (40) and (41).

Since $|\tau_{j}|\leq 1$ and using (42), we have a.s. that for $n>n(\omega)$,

The norm convergence (10) implies that $\sum_{j}Y_{n}(j)\rightarrow 2c_{0}\sum\varrho_{j}$ and so it follows from the version of the dominated convergence theorem due to Pratt (1960) that

Proof of Theorem 3 [Theorem 1 is a special case.] We apply the perturbation theorem with $A=R$ and $E=A+B+C$. The separation between the principal eigenvalue of $R$ and the remaining ones is

while from Proposition 1 we have the bound

A.3 Lower Bounds: Proof of Theorem 2

We begin with a heuristic outline of the proof. We write $S$ in the form $D+B$, introducing

while, as before, $B=\rho u^{T}+u\rho^{T}$ and $u=\sigma n^{-1}Zv$.

A symmetry trick plays a major role: write $S_{-}=D-B$ and let $\hat{\rho}_{-}$ be the principal unit eigenvector for $S_{-}$.

The argument makes precise the following chain of remarks, which are made plausible by reference to Figure 8.

(i) $B\hat{\rho}$ is nearly orthogonal to $D\hat{\rho}+B\hat{\rho}=S\hat{\rho}=\hat{\lambda}\hat{\rho}$.

(ii) the side length $\|B\hat{\rho}\|$ is bounded away from zero, when $c>0$.

(iii) the angle between $\hat{\rho}$ and $S_{-}\hat{\rho}$ is “large”, i.e. bounded away from zero.

(iv) the angle between $\hat{\rho}$ and $\hat{\rho}_{-}$ is “large” [this follows from Lemma 1 applied to $M=S_{-}$.]

(v) and finally, the angle between $\hat{\rho}$ and $\rho$ must be “large”, due to the equality in distribution of $\hat{\rho}$ and $\hat{\rho}_{-}$.

Getting down to details, we will establish (i)-(iii) under the assumption that $\hat{\rho}$ is close to $\rho$. Specifically, we show that given $\delta>0$ small, there exists $\alpha(\delta)=\alpha(\delta;\sigma,c)>0$ such that w.p. $\rightarrow 1$,

(ii’) $\|Bx\|$ is bounded below (see (49)), and

(i’) $Bx$ is nearly orthogonal to $x$ (see (50).

For convenience in this proof, we may take $\|\rho\|=1$. Write $x\in N_{\delta_{1}}$ in the form

Since $B\rho=(u^{T}\rho)\rho+u$ and $B\eta=(u^{T}\eta)\rho$, we find that

Denote the second right side term by $r$: clearly $\|r\|\leq|u^{T}\rho|+(\sin\delta)\|u\|$, and so, uniformly on $N_{\delta}$,

Since both $\|u\|\rightarrow\sigma\sqrt{c}$ and $u^{T}\rho\rightarrow 0$ a.s., we conclude that w.p. $\rightarrow 1$,

Turning to the angle between $x$ and $Bx$, we find from (47) and (48) that

Consequently, using $\|x\|=1$ and (49), w.p. $\rightarrow 1$, and for $\delta<\pi/4$, say,

Now return to Figure 8. As a prelude to step (iii), we establish a lower bound for $\alpha=\angle(\hat{\rho},D\hat{\rho}).$ Applying the sine rule (23) to $\xi=D\hat{\rho}$ and $\eta=\hat{\lambda}\hat{\rho}=D\hat{\rho}+B\hat{\rho}$, we obtain

On the assumption that $\hat{\rho}\in N_{\delta}$, bound (50) yields

On the other hand, since $\|\hat{\rho}\|=1$,

Combining the last three bounds into (51) shows that there exists a positive $\alpha(\delta;\sigma,c)$ such that if $\hat{\rho}\in N_{\delta}$, then w.p. 1 for large $n$,

Returning to Figure 8, consider $\angle(D\hat{\rho}+B\hat{\rho},D\hat{\rho}-B\hat{\rho})=\alpha+\gamma$. Since $\beta\geq\pi/2-c_{3}\delta$, we clearly have $\alpha+\gamma\leq\pi-\beta\leq\pi/2+c_{3}\delta$ and hence

In particular, with $\delta\leq\delta_{0}(\sigma,c)$,

which is our step (iii). As mentioned earlier, Lemma 1 applied to $M=S_{-}$ entails that $\angle(\hat{\rho},\hat{\rho}_{-})\geq(1/3)\alpha(\delta)$. For the rest of the proof, we write $\hat{\rho}_{+}$ for $\hat{\rho}$. To summarize to this point, we have shown that if $\angle(\hat{\rho}_{+},\rho)\leq\delta$, then w.p. $\rightarrow 1$,

Note that $S$ and $S_{-}$ have the same distribution: viewed as functions of random terms $Z$ and $v$:

We call an event $\mathcal{A}$ symmetric if $(Z,v)\in\mathcal{A}$ iff $(Z,-v)\in\mathcal{A}$. For such symmetric events

From this and the triangle inequality for angles

By the symmetry of the distributions, conclusion (52) is also obtained w.p. $\rightarrow 1$ if $\angle(\hat{\rho}_{-},\rho)\leq\delta$. Consequently, letting $\mathcal{A}$ refer to the symmetric event $\mathcal{A}_{\delta}=\{\angle(\hat{\rho}_{+},\rho)\leq\delta\}\cup\{\angle(\hat{\rho}_{-},\rho)\leq\delta\}$, we have

This completes the proof of Theorem 2. The lower bound proof for Theorem 3 proceeds similarly, but is omitted – for some extra detail, see Lu (2002).

A.4 Proof of Theorem 4.

We may assume, without loss of generality, that $\sigma_{1}^{2}\geq\sigma_{2}^{2}\geq\cdots\geq\sigma_{p}^{2}.$

False inclusion. For any fixed constant $t$,

This threshold device leads to bounds on error probabilities using only marginal distributions. For example, consider false inclusion of variable $l$:

Write $\bar{M}_{n}$ for a $\chi^{2}_{(n)}/n$ variate, and note from (15) that $\hat{\sigma}_{\nu}^{2}\sim\sigma_{\nu}^{2}\bar{M}_{n}$. Set $t=\sigma_{k}^{2}(1-\epsilon_{n})$ for a value of $\epsilon_{n}$ to be determined. Since $\sigma_{i}^{2}\geq\sigma_{k}^{2}$ and $\sigma_{l}^{2}\leq\sigma_{k}^{2}(1-\alpha_{n})$, we arrive at

using large deviation bound (29). With the choice $\epsilon_{n}=\sqrt{3}\alpha_{n}/(2+\sqrt{3})$, both exponents are bounded above by $-b(\gamma)\log n$, and so $P\{FI\}\leq p(k+1)n^{-b(\gamma)}$.

False exclusion. The argument is similar, starting with the remark that for any fixed $t$,

Consequently, if we set $t=\sigma_{k}^{2}(1+\epsilon_{n})$ and use $\sigma_{l}^{2}\geq\sigma_{k}^{2}(1+\alpha_{n})$, we get

The bound $P\{FE\}\leq pkn^{-b(\gamma)}+ke^{-b(\gamma)(1-2\alpha_{n})\log n}$ follows on setting $\epsilon_{n}=2\alpha_{n}/(2+\sqrt{3})$ and noting that $(1+\alpha_{n})^{-2}\geq 1-2\alpha_{n}$ .

For numerical bounds, we may collect the preceding bounds in the form

A.5 Proof of Theorem 5

Outline. Recall that $\gamma_{n}=\gamma(n^{-1}\log n)^{1/2},$ that the selected subset of variables $\hat{I}$ is defined by

and that the estimated principal eigenvector based on $\hat{I}$ is written $\hat{\rho}_{I}$. We set

and will use the triangle inequality $d(\hat{\rho}_{I},\rho)\leq d(\hat{\rho}_{I},\rho_{I})+d(\rho_{I},\rho)$ to show that $\hat{\rho}_{I}\rightarrow\rho$. There are three main steps.

(i) Construct deterministic sets of indices

which bracket $\hat{I}$ almost surely as $n\rightarrow\infty$:

(ii) the uniform sparsity, combined with $\hat{I}^{c}\subset I_{n}^{-c},$ is used to show that

(iii) the containment $\hat{I}\subset I_{n}^{+}$, combined with $|I_{n}^{+}|=o(n)$ shows via methods similar to Theorem 4 that

Details. Step (i). We first obtain a bound on the cardinality of $I_{n}^{\pm}$ using the uniform sparsity conditions (17). Since $|\rho|_{(\nu)}\leq C\nu^{-1/q}$

Turning to the bracketing relations (55), we first remark that $\hat{\sigma}_{\nu}^{2}\stackrel{{\scriptstyle\mathcal{D}}}{{=}}\sigma_{\nu}^{2}\chi_{(n)}^{2}/n$, and when $\nu\in I_{n}^{\pm}$,

Using the definitions of $\hat{I}$ and writing $\bar{M}_{n}$ for a random variable with the distribution of $\chi^{2}_{(n)}/n$, we have

We apply (29) with $\epsilon_{n}=(a_{+}-1)\gamma_{n}/(1+a_{+}\gamma_{n})$ and for $n$ large and $\gamma^{\prime}$ slightly smaller than $\gamma^{2}$,

with $\gamma_{+}^{{}^{\prime\prime}}=(a_{+}-1)^{2}\gamma^{\prime}/4.$ If $\sqrt{\gamma}\geq 12,$ then $\gamma_{+}^{{}^{\prime\prime}}\geq 3$ for suitable $a_{+}>2.$

The argument for the other inclusion is analogous:

with $\gamma_{-}^{{}^{\prime\prime}}=3(1-a_{-})^{2}\gamma^{\prime}/16$ so long as $n$ is large enough. If $\sqrt{\gamma}\geq 12$, then $\gamma_{-}^{{}^{\prime\prime}}>2$ for suitable $a_{-}<1-\sqrt{8/9}$.

By a Borel-Cantelli argument, (55) follows from the bounds on $P_{n}^{-}$ and $P_{n}^{+}$.

Step (ii). For $n>n(\omega)$ we have $I_{n}^{-}\subset\hat{I}$ and so

When $\nu\in I_{n}^{-c}$, we have by definition

say, while the uniform sparsity condition entails

Putting these together, and defining $s_{*}=s_{*}(n)$ as the solution of the equation $Cs^{-1/q}=\epsilon_{n}$, we obtain

As in the proof of Theorem 3, we consider $S^{*}_{I}=S_{I}-\sigma^{2}I_{\hat{k}}=\rho_{I}\rho^{T}_{I}+E_{I}$ and note that the perturbation term has the decomposition

Consider the first term on the right side. Since $\|\rho_{I}-\rho\|_{2}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}0$ from step (ii), it follows that $\|\rho_{I}\|_{2}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}\|\rho\|$. As before $v_{s}\stackrel{{\scriptstyle a.s.}}{{\rightarrow}}0$, and so the first term is asymptotically negligible.

Let $Z_{I^{+}}=(z_{\nu i}~{}:~{}\nu\in I^{+}_{n},i=1,\ldots,n)$ and $u_{I^{+}}=(u_{\nu}~{}:~{}\nu\in I_{n}^{+})$. On the event $\Omega_{n}=\{\hat{I}\subset I_{n}^{+}\}$, we have

and setting $k_{+}=|I_{n}^{+}|$, by the same arguments as led to (40), we have

Finally, since on the event $\Omega_{n}$, the matrix $Z_{I^{+}}$ contains $Z_{I}$, along with some additional rows, it follows that

by (38), again since $k_{+}=o(n)$. Combining the previous bounds, we conclude that $\|E_{I}\|_{2}\rightarrow 0$.

The separation $\delta_{n}=\|\rho_{I}\|_{2}^{2}\rightarrow\|\rho\|_{2}^{2}>0$ and so by the perturbation bound

Acknowledgements. The authors are grateful for helpful comments from Debashis Paul and the participants at the Functional Data Analysis meeting at Gainesville, FL. January 9-11, 2003. This work was supported in part by grants NSF DMS 0072661 and NIH EB R01 EB001988.

References