Augmented sparse principal component analysis for high dimensional data

Introduction

Principal components analysis (PCA) has been a widely used technique in reducing dimensionality of multivariate data. A traditional setting where PCA is applicable is when one has repeated observations from a multivariate population that can be described reasonably well by its first two moments. When the dimension of sample observations, is fixed, distributional properties of the eigenvalues and eigenvectors of the sample covariance have been dealt with at length by various authors. Anderson (1963), Muirhead (1982) and Tyler (1983) are among standard references. Much of the “large sample” study of the eigen-structure of the sample covariance matrix is based on the fact that, sample covariance approximates population covariance matrix well when sample size is large. However, due to advances in data acquisition technologies, statistical problems, where the dimensionality of individuals are of nearly the same order of magnitude as (or even bigger than) the sample size, are increasingly common. The following is a representative list of areas and articles where PCA has been in use. In all these cases $N$ denotes the dimension of an observation and $n$ denotes the sample size.

Image recognition : The face recognition problem is to identify a face from a collection of faces. Here each observation is a digitized image of the face of a person. So typically, with $128\times 128$ pixel grids, one has to deal with a situation where $N\approx 1.6\times 10^{6}$. Whereas, a standard image database, e.g. that of students of Brown University Wickerhauser (1994), may contain only a few hundred pictures.

Shape analysis : Stegmann and Gomez (2002), Cootes, Edwards and Taylor (2001) outline a class of methods for analyzing the shape of an object based on repeated measurements that involves annotating the objects for landmarks. These landmarks act as features of the objects, and hence, can be thought of as the dimension of the observations. For a specific example relating to motion of hand Stegmann and Gomez (2002), the number of landmarks is 56 and sample size is 40.

Chemometrics : In many chemometric studies, sometimes the data consists of several thousands of spectra measured at several hundred wavelength positions, e.g. data collected for calibration of spectrometers. Vogt, Dable, Cramer and Booksh (2004) give an overview of some of these applications.

Econometrics : Large factor analysis models are often used in econometric studies, e.g. in dealing with hundreds of stock prices as a multivariate time series. Markowitz’s theory of optimal portfolios ask this question. Given a set of financial assets characterized by their average return and risk, what is the optimal weight of each asset, such that the overall portfolio provides the best return? Laloux, Cizeau, Bouchaud and Potters (2000) discuss several applications. Bai (2003) considers some inferential aspects.

Climate studies : Measurements on atmospheric indicators, like ozone concentration etc. are taken at a number of monitoring stations over a number of time points. In this literature, principal components are commonly referred to as “empirical orthogonal functions”. Preisendorfer (1988) gives a detailed treatment. EOFs are also used for model diagnostics and data summary Cassou, Deser, Terraty, Hurrell and Drévillon (2004).

Communication theory : Tulino and Verdu (2004) give an extensive treatment to the connection between random matrix theory and vector channels used in wireless communications.

Functional data analysis : Since observations are curves, which are typically measured at a large number of points, the data is high dimensional. Buja, Hastie and Tibshirani (1995) give an example of speech dataset consisting of 162 observations - each one is a periodogram of a “phoneme” spoken by a person. Ramsay and Silverman (2002) discuss other applications.

Microarray analysis : Gene microarrays present data in the form expression profiles of several thousand genes for each subject under study. Bair, Hastie, Paul and Tibshirani (2006) analyze an example involving the study of survival times of 240 ($=n$) patients with diffuse large B-cell lymphoma, with gene expression measurements for 7389 ($=N$) genes.

Of late, researchers in various fields have been using different versions of non-identity covariance matrices of growing dimension. Among these, a particularly interesting model assumes that,

the eigenvalues of the population covariance matrix $\Sigma$ are (in descending order)

This has been deemed the “spiked population model” by Johnstone (2001). It has also been observed that for certain types of data, e.g. in speech recognition Buja, Hastie and Tibshirani (1995), wireless communication Telatar (1999), statistical learning (Hoyle and Rattray (2003, 2004)), a few of the sample eigenvalues have limiting behavior that is different from the behavior when the covariance is the identity. This paper deals with the issue of estimating the eigenvectors of $\Sigma$, when it has the structure described by (*), and the dimension $N$ grows to infinity together with sample size $n$.

In many practical problems, at least the leading eigenvectors are thought to represent some underlying phenomena. This has been one of the reasons for their popularity in analysis of what can be characterized as functional data. For example, Zhao, Marron and Wells (2004) consider the “yeast cell cycle” data of Spellman et al. (1998), and argue that the first two components obtained by a functional PCA of the data represent systematic structure. In climate studies, empirical orthogonal functions are often used for identifying patterns in the data, as well as for data summary. See for example Corti, Molteni and Palmer (1999). In many of these instances there is some idea about the structure of the eigenvectors of the covariance matrix, such as to the extent they are smooth, or oscillatory. At the same time, these data are often corrupted with a substantial amount of noise, which can lead to very noisy estimates of the eigen-elements. There is also a growing literature on functional response models in which the regressors are random functions and the responses are either vectors or functions (Chiou, Müller and Wang (2004), Hall and Horowitz (2004), Cardot, Ferraty and Sarda (2003)). Quite often a functional principal component regression is used to solve these problems. Thus, there are both practical and scientific interests in devising methods for estimating the eigenvectors and eigenvalues that can take advantage of the information about the structure of the population eigenvectors. At the same time, there is also a need to address this estimation problem from a broader statistical perspective.

In multivariate analysis, there is a huge body of work on estimation of population covariance, and in particular on developing optimal strategies for estimation from a decision theoretic point of view. Dey and Srinivasan (1985), Efron and Morris (1976), Haff (1980), Loh (1988) are some of the standard references in this field. However, a decision theoretic treatment of functional data analysis is still somewhat limited in its breadth. Hall and Horowitz (2004) and Tony Cai and Hall (2005) derive optimal rates of convergence of estimators of the regression function and fitted response in functional linear model context. Cardot (2000) gave upper bounds on the rate of convergence of a spline-based estimator of eigenvectors under some smoothness assumptions. Kneip (1994) also derived similar results in a slightly different context.

In this paper, the aim is to address the problem of estimating eigenvectors from a minimax risk analysis viewpoint. Henceforth, the observations will be assumed to have a Gaussian distribution. This assumption, though somewhat idealized, helps in bringing out some essential features of the estimation problem. Since algebraic manipulation of spectral elements of a matrix is rather difficult, it is not easy to make any precise finite sample statement about the risk properties of estimators. Therefore the analysis is mostly asymptotic in nature, even though efforts have been made to make the approximations to risk etc. as explicit as possible. The asymptotic regime considered here assumes a triangular array structure in which $N$, the dimensionality of individual observations, tends to $\infty$ with sample size $n$. This framework is partly motivated by similar analytical approaches to the problem of estimation of mean function in nonparametric regression context. In particular, a squared error type loss is proposed, and some $l^{q}$-type sparsity constraint is imposed on the parameters, which in our case are individual eigenvectors. Relevance of this sort of constraints in the context of functional data analysis is discussed in Section 3. The main results of this chapter are the following. Theorem 1 describes risk behavior of sample eigenvectors as estimators of their population counterparts. Theorem 2 gives a lower bound on the minimax risk. An estimation scheme, named Augmented Sparse Principal Component Analysis (ASPCA) is proposed and is shown to have the optimal rate of convergence over a class of $l^{q}$ norm-constrained parameter spaces under suitable regularity conditions. Throughout it is assumed that the leading eigenvalues of the population covariance matrix are distinct, so the eigenvectors are identifiable. A more general framework, which looks at estimating the eigen-subspaces and allows for eigenvalues with arbitrary multiplicity, is beyond the scope of this paper.

Model

Suppose that, $\{X_{i}:i=1,\ldots,n\}_{n\geq 1}$ is a triangular array, where the $N\times 1$ vectors $X_{i}:=X_{i}^{n},i=1,\ldots,n$ are i.i.d. on a common probability space for each $n$. The dimension $N$ is assumed to be a function of $n$ and increases without bound as $n\to\infty$. The observation vectors are assumed to be i.i.d. as $N(\xi,\Sigma)$, where $\xi$ is the mean vector; and $\Sigma$ is the covariance matrix. The assumption on $\Sigma$ is that, it is a finite rank perturbation of (a multiple of) the identity. In other words,

where $\lambda_{1}>\lambda_{2}>\ldots>\lambda_{M}>0$, and the vectors $\theta_{1},\ldots,\theta_{M}$ are orthonormal. Notice that strict inequality in the order relationship among the $\lambda_{\nu}$’s implies that the $\theta_{\nu}$ are identifiable up to a sign convention. Notice that with this identifiability condition, $\theta_{\nu}$ is the eigenvector corresponding to the $\nu$-th largest eigenvalue, namely, $\lambda_{\nu}+\sigma^{2}$, of $\Sigma$. The term “finite rank” means that, $M$ will remain fixed for all the asymptotic analysis that follows. This analysis involves letting both $n$ and $N$ increase to infinity simultaneously. Therefore, $\Sigma$, the $\lambda_{\nu}$’s and the $\theta_{\nu}$’s should be thought of as being dependent on $N$.

The observations can be equivalently described in terms of the factor analysis model :

Here, for each $n$, $v_{\nu i}$, $Z_{ik}$ are all independently and identically distributed as $N(0,1)$. $M\geq 1$ is assumed fixed.

Since the eigenvectors of $\Sigma$ are invariant to a scale change in the original observations, for simplifying notation, it is assumed that $\sigma=1$. Notice that this also means that, $\lambda_{1},\ldots,\lambda_{M}$ appearing in the results relating to the rates of convergence of various estimators of $\theta_{\nu}$ should be changed to $\lambda_{1}/\sigma,\ldots,\lambda_{M}/\sigma$ when (1) holds with an arbitrary $\sigma>0$.

Another simplifying assumption is that, $\xi=0$. This is because, the main focus of the current exposition is on estimating the eigen-structure of $\Sigma$, and the unnormalized sample covariance matrix

where $\overline{X}$ is the sample mean, has the same distribution as that of the matrix

where $Y_{i}$ are i.i.d. $N(0,\Sigma)$. This means that, for estimation purposes, if the attention is restricted to the sample covariance matrix, then from an asymptotic analysis point of view, it is enough to assume $\xi=0$, and to define the sample covariance matrix as $\mathbf{S}=\frac{1}{n}\mathbf{X}\mathbf{X}^{T}$, where $\mathbf{X}=[X_{1}:\ldots:X_{n}]$.

The following condition, or Basic Assumption will be used frequently, and will be referred to as BA.

(2) and (1) hold, with $\xi=0$ and $\sigma=1$; $N=N(n)\to\infty$ as $n\to\infty$; $\lambda_{1}>\ldots>\lambda_{M}>0$.

The goal is, given data $X_{1},X_{2},\ldots,X_{n}$, to estimate $\theta_{\nu}$, for $\nu=1,\ldots,M$. To assess the performance of any such estimator, a minimax risk analysis approach is proposed. The first task is to specify a loss function for this estimation problem. Observe that since the model is invariant under separate changes of sign of the $\theta_{\nu}$, it is necessary to specify a loss function that is also invariant under a sign change. We specify the following loss function :

2 Rate of convergence for ordinary PCA

It is assumed that either $\lambda_{1}$ is fixed, or that it varies with $n$ and $N$ so that,

as $n,N\to\infty$, $\frac{\lambda_{\nu}}{\lambda_{1}}\to\rho_{\nu}$ for $\nu=1,\ldots,M$, where $1=\rho_{1}>\rho_{2}>\ldots>\rho_{M}$;

as $n,N\to\infty$, $\frac{N}{nh(\lambda_{1})}\to 0$, where

Notice that, all four conditions (i)-(iv) below imply that $\frac{N}{nh(\lambda_{1})}\to 0$ as $n\to\infty$.

$\frac{N}{n}\to\gamma\in(0,\infty)$ and $\frac{N}{n\lambda_{1}}\to 0$

$\lambda_{1}\to 0$, $\frac{N}{n}\to 0$ and $\frac{N}{n\lambda_{1}^{2}}\to 0$

$0<\lim\inf_{n\to\infty}\lambda_{1}\leq\lim\sup_{n\to\infty}\lambda_{1}<\infty$ and $\frac{N}{n}\to 0$

$\frac{N}{n}\to\infty$, and $\frac{N}{n\lambda_{1}}\to 0$.

Remark : Condition L1 is really an asymptotic identifiability condition which guarantees that at the scale of the largest “signal” eigenvalue, bigger eigenvalues are well-separated.

Theorem 1: Suppose that the eigenvalues $\lambda_{1},\ldots,\lambda_{M}$ satisfy L1 and L2. If $\log(n\vee N)=o(n\wedge N)$, then for $\nu=1,\ldots,M$,

Remark : It is possible to relax some of the conditions stated in the theorem. On the other hand, with some reasonable assumptions on the decay of the eigenvalues, it is also possible to incorporate cases where $M$ is no longer a constant, but increases with $n$. Then the issues would include, rates of growth of $M$ and the rate of decay of eigenvalues that would result in the OPCA estimator retaining consistency and the expression for its asymptotic risk. These issues are not going to be addressed here. However, it is important to note that, such questions have been investigated - not necessarily for the Gaussian case - in the context of spectral decomposition of $L^{2}$ stochastic processes by, Hall and Horowitz (2004), Tony Cai and Hall (2005), Boente and Fraiman (2000), Hall and Hosseini-Nasab (2006) among others. However, these analyses do not deal with measurement errors. The condition $\frac{N}{nh(\lambda_{\nu})}\to 0$ is a necessary condition for uniform convergence, as shown in Theorem 2. It should be noted that, there are results, proved under slightly different circumstances, that obtain the rates given by (5) as an upper bound on the rate of convergence of OPCA estimators (Bai (2003), Cardot (2000), Kneip (1994)). These analyses, while treating the problem under less restrictive assumptions than Gaussianity (essentially, finite eighth moment for the noise $Z_{ik}$), make the assumption that $\frac{N^{2}}{n}\to 0$, when the $\lambda_{\nu}$’s are considered fixed.

Sparse model for eigenvectors

In this section we discuss the concept of sparsity of the eigenvectors and impose some restrictions on the space of eigenvectors that lead to a sparse parametrization. This notion will be used later from a decision-theoretic view point in order to analyze the risk behavior of estimators of the eigenvectors. From now on, $\theta$ will be used to denote the matrix $[\theta_{1},\ldots,\theta_{M}]$.

The parameter space is taken to be a class of $M$-dimensional positive semi-definite matrices satisfying the following criteria:

$\theta_{1},\ldots,\theta_{M}$ are orthonormal.

In the Functional Data Analysis context, one can think of the observations as the vectors of wavelet coefficients (when transformed in an orthogonal wavelet basis of sufficient regularity) of the observed functions. If the smoothness of a function $g$ is measured by its membership in a Besov space $B_{q^{\prime},r}^{\alpha}$, and if the vector of its wavelet coefficients, when expanded in a sufficiently regular wavelet basis, is denoted by $\mathbf{g}$, then from Donoho (1993),

One may refer to Johnstone (2002) for more details. Treating this as a motivation, instead of imposing a weak-$l^{q}$ constraint on the parameter $\theta_{\nu}$, we rather impose an $l^{q}$ constraint. Note that, for $C,q>0$,

Since $l^{q}(C)\hookrightarrow wl^{q}(C)$, it is possible to derive lower bounds on the minimax risk of estimators when the parameter lies in a $wl^{q}$ space by restricting attention to an $l^{q}$ ball of appropriate radius.

Then $m_{C}$ is the largest dimension of a unit sphere, centered at 0, that fits inside the parameter space $\Theta_{q}(C)$.

2 Parameter space

The parameter space for $\theta:=[\theta_{1}:\ldots:\theta_{M}]$ is denoted by

where $\Theta_{q}(C)$ is defined through (8), and $C_{\nu}\geq 1$ for all $\nu=1,\ldots,M$.

Remark : If $M>1$, one can describe the sparsity of the eigenvectors in a different way. Consider the sequence $\zeta:=\zeta^{N}=(\sqrt{\sum_{\nu=1}^{M}\lambda_{\nu}\theta_{\nu k}^{2}}:k=1,2,\ldots,N)$. One may demand that the vector $\zeta$ be sparse in an $l^{q}$ or weak-$l^{q}$ sense. This particular approach to sparsity has some natural interpretability, since the quantity $\zeta_{k}^{2}=\sum_{\nu=1}^{M}\lambda_{\nu}\theta_{\nu k}^{2}$, where $\zeta_{k}$ is the $k$-th coordinate of $\zeta$, is the variance of the $k$-th coordinate of the “signal” part of the vector $X$. There is a connection between this model and the model we intend to study. If (10) holds, then $\zeta\in l^{q}_{N}(\overline{C}_{\lambda})$, where $\overline{C}_{\lambda}^{q}=\sum_{\nu=1}^{M}\lambda_{\nu}^{q/2}C_{\nu}^{q}$. On the other hand, $l^{q}$ (weak-$l^{q}$) sparsity of $\zeta$ implies $l^{q}$ (weak-$l^{q}$) sparsity of $\theta_{\nu}$ for all $\nu=1,\ldots,M$.

3 Lower bound on the minimax risk

In this section a lower bound on the minimax risk of estimating $\theta_{\nu}$ over the parameter space (10) is derived when $0<q<2$, under the loss function defined through (3). The result is stated under some simplifying assumptions that make the asymptotic analysis more transparent. Define

There exists a constant $C_{0}>0$ such that $C_{0}^{q}<C_{\mu}^{q}-1$ for all $\mu=1,\ldots,M$, for all $N$.

As $n,N\to\infty$, $nh(\lambda_{\nu})\to\infty$.

As $n,N\to\infty$, $nh(\lambda_{\nu})=O(1)$.

As $n,N\to\infty$, $ng(\lambda_{\mu},\lambda_{\nu})\to\infty$ for all $\mu=1,\ldots,\nu-1,\nu+1,\ldots,M$.

As $n,N\to\infty$, $n\max_{1\leq\mu\neq\nu\leq M}g(\lambda_{\mu},\lambda_{\nu})=O(1)$.

Conditions A4 and A5 are applicable only when $M>1$. In the statement of the following theorem, the infimum is taken over all estimators $\widehat{\theta}_{\nu}$, estimating $\theta_{\nu}$, satisfying $\parallel\widehat{\theta}_{\nu}\parallel=1$.

Theorem 2: Let $0<q<2$ and $1\leq\nu\leq M$. Suppose that A1 holds.

If A3 holds, then there exists $B_{1}>0$ such that

If A2 holds, then there exists $B_{2}>0$, $A_{q}>0$, and $c_{1}\in(0,1)$, such that

for some $K>0$, $\alpha\in(0,1)$, $c_{q}(\alpha)\in(0,1)$ and $A_{q,\alpha}>0$. Here $\overline{C}_{\nu}^{q}:=C_{\nu}^{q}-1$. Also, one can take $c_{1}=\log(9/8)$, $A_{q}=(\frac{9c_{1}}{2})^{1-q/2}$, $A_{q,\alpha}=(\alpha/2)^{1-q/2}$, $c_{q}(\alpha)=(\alpha/9)^{1-q/2}$, $B_{2}=\frac{1}{8}$ and $B_{3}=(8e)^{-1}$.

Suppose that $M>1$. If A4 holds, then there exists $B_{3}>0$ such that

One can take $B_{3}=\frac{1}{8e}$. However, if A5 holds, then (12) is true.

Remark : In the statement of Theorem 2, there is much flexibility in terms of what values the “hyperparameters” $C_{1},\ldots,C_{M}$ and the eigenvalues $\lambda_{1},\ldots,\lambda_{M}$ can take. In particular, they can vary with $N$, subject to the modest requirement that A1 is satisfied. However, the constants appearing in equations (13) and (16) are not optimal.

Remark : Another notable aspect is that, as the proof later shows, the rate lower bounds in Part (b) are all of the form $\frac{m}{nh(\lambda_{\nu})}$, where $m$ is the “effective” number of “significant” coordinates. This phrase becomes clear if one notices further that, in the construction that leads to the lower bound (see Section 6.7), the vector $\theta_{\nu}$ in a near-worst case scenario has overwhelming number of coordinates of size $const.~\frac{1}{\sqrt{nh(\lambda_{\nu})}}$, or, in the case (15), of size $const.~\frac{\sqrt{\log N}}{\sqrt{nh(\lambda_{\nu})}}$. Here $m$ is of the same order as the number of these “significant” coordinates. This suggests that, an estimation strategy that is able to extract coordinates of $\theta_{\nu}$ of the stated size, would have the right rate of convergence, subject to possibly some regularity conditions. The estimator described later (ASPCA) is constructed by following this principle.

Part (a) and the second statement of Part (c) of Theorem 2 depict situations under which there is no estimator that is asymptotically uniformly consistent over $\Theta_{q}^{M}(C_{1},\ldots,C_{M})$. Moreover, the first part of Part (b), and Theorem 1 readily yield the following corollary.

Corollary 1: If the conditions of Theorem 1 hold, and if A1 holds, together with the condition that

then the usual PCA-based estimator of $\widehat{\theta}_{\nu}$, i.e. the eigenvector corresponding to the $\nu$-th largest eigenvalue of $\mathbf{S}$, has asymptotically the best rate of convergence.

Remark : A closer look at the proof of Theorem 1 reveals that the method of proof explicitly made use of condition L1 to ensure that the contribution of $\lambda_{1},\ldots,\lambda_{M}$ to the residual term of the second order expansion of $\widehat{\theta}_{\nu}$ is bounded. However, the condition $n\max_{\mu\neq\nu}g(\lambda_{\mu},\lambda_{\nu})\to\infty$ is certainly much weaker than that. The method of proof pursued here fails to settle the question as to whether this is sufficient to get the asymptotic rate (5). It is conjectured that this is the case.

Estimation scheme

This section outlines an estimation strategy for the eigenvectors $\theta_{\nu}$, $\nu=1,\ldots,M$. Model (2) is assumed throughtout for observations $X_{i}$, $i=1,\ldots,n$. We propose estimators is for the case when the noise variance $\sigma^{2}$ is known. Therefore, without loss of generality, it can be taken to be 1. Henceforth, for simplicity of notations, it is also assumed that $\xi=0$. In practice, one may have to estimate $\sigma^{2}$ from data. The median of the diagonal entries of the sample covariance matrix $\mathbf{S}:=\frac{1}{n}\mathbf{X}\mathbf{X}^{T}$ serves as a reasonable (although slightly biased) estimator of $\sigma^{2}$, if the true model is sparse. In the latter case, the data are rescaled by multiplying each observation by $\widehat{\sigma}^{-1}$, and the resultant covariance matrix is called, with a slight abuse of notation, $\mathbf{S}$. Note that, in this case, the estimates of eigenvalues of $\Sigma$ are $\widehat{\sigma}^{2}$ times the corresponding eigenvalues of $\mathbf{S}$.

In order to motivate the approach that is described in what follows, consider first the SPCA estimation scheme studied by Johnstone and Lu (2004). To that end, let $\mathbf{S}=\frac{1}{n}\mathbf{X}\mathbf{X}^{T}$ denote the sample covariance matrix. Suppose that the sample variances of coordinates (i.e., diagonal terms of $\mathbf{S}$) are denoted by $\hat{\sigma}_{1}^{2},\ldots,\hat{\sigma}_{N}^{2}$.

Define $\widehat{I}_{n}$ to be the set of indices $k\in\{1,\ldots,N\}$ such that $\hat{\sigma}_{k}^{2}>\gamma_{n}$ for some threshold $\gamma_{n}>0$.

Let $\mathbf{S}_{\widehat{I}_{n},\widehat{I}_{n}}$ be the submatrix of $\mathbf{S}$ corresponding to the coordinates $\widehat{I}_{n}$. Perform an eigen-analysis of $\mathbf{S}_{\widehat{I}_{n},\widehat{I}_{n}}$. Denote the eigenvectors by $\mathbf{e}_{1},\ldots,\mathbf{e}_{\min\{n,|\widehat{I}_{n}|\}}$.

For $\nu=1,\ldots,M$, estimate $\theta_{\nu}$ by $\widetilde{\mathbf{e}}_{\nu}$ where $\widetilde{\mathbf{e}}_{\nu}$, an $N\times 1$ vector, is obtained from $\mathbf{e}_{\nu}$ by augmenting zeros to all the coordinates that are in $\{1,\ldots,N\}\setminus\widehat{I}_{n}$.

Johnstone and Lu (2004) showed that, if one chooses an appropriate threshold $\gamma_{n}$, then the estimate of $\theta_{\nu}$ is consistent under the weak-$l^{q}$ sparsity constraint on $\theta_{\nu}$. However, Paul and Johnstone (2004) showed that even with the best choice of $\gamma_{n}$, the rate of convergence of the risk of this estimate is not optimal. Indeed, Paul and Johnstone (2004) demonstrate an estimator which has a better rate of convergence in the single component ($M=1$) situation.

2 Augmented Sparse PCA (ASPCA)

We now propose the ASPCA estimation scheme. This scheme is a refinement of the SPCA scheme of Johnstone and Lu (2004), and can be viewed as a generalization of the estimation scheme proposed by Paul and Johnstone (2004) in the single component ($M=1$) case.

The key idea behind this estimation scheme is that, in addition to using the coordinates having large variance, if one also uses the covariance structure appropriately, then under the assumption of a sparse structure of the eigenvectors, one will be able to extract a lot more information and thereby get more accurate estimate of the eigenvalues and eigenvectors. Notice that SPCA only focuses on the diagonal of the covariance matrix and therefore ignores the covariance structure. This renders this scheme suboptimal from an asymptotic minimax risk analysis point of view. To make this point clearer, it is instructive to analyze the covariance matrix in the $M=1$ case. In view of the second Remark after the statement of Theorem 2 one expects to be able to recover coordinates $k$ for which $|\theta_{1k}|\gg\frac{1}{\sqrt{nh(\lambda_{1})}}$. However, the best choice for $\gamma_{n}$ for SPCA is $\gamma\sqrt{\frac{\log n}{n}}$, for some constant $\gamma>0$, which is way too large. On the other hand, suppose that one divides the coordinates into two sets $A$ and $B$, where the former contains all those $k$ such that $|\theta_{k}|$ is “large”, and the latter contains smaller coordinates. Partition the matrix $\Sigma$ as

Here $\Sigma_{BA}=\lambda_{1}\theta_{1,B}\theta_{1,A}^{T}$. Assume that, there is a “preliminary” estimator of $\theta_{1}$, say $\widetilde{\theta}_{1}$ such that, $\langle\widetilde{\theta}_{1,A},\theta_{1,A}\rangle\to 1$ in probability as $n\to\infty$. Then one can use this estimator as a “filter”, in a way described below, to recover the “informative ones” among the smaller coordinates. This can be seen from the following relationship

In this manner one can extract some information about those coordinates of $\theta_{1}$ that are in set $B$. The algorithm described below is a generalization of this idea. It has three stages. First two stages will be referred to as “coordinate selection” stages. The final stage consists of an eigen-analysis of the submatrix of $\mathbf{S}$ corresponding to the selected coordinates, followed by a hard thresholding of the estimated eigenvectors.

Let $\gamma_{i}>0$ for $i=1,2,3$ and $\kappa>0$ be four constants to be specified later. Define $\gamma_{1,n}=\gamma_{1}\sqrt{\frac{\log(n\vee N)}{n}}$.

Select coordinates $k$ such that $\widehat{\sigma}_{kk}:=\mathbf{S}_{kk}>1+\gamma_{1,n}$. Denote the set of selected coordinates by $\widehat{I}_{1,n}$.

Denote the diagonal of the matrix $\mathbf{Q}\mathbf{Q}^{T}$ by $T$. Define $\widehat{I}_{2,n}$ to be the set of coordinates $k\in\{1,\ldots,N\}\setminus\widehat{I}_{1,n}$ such that $|T_{k}|>\gamma_{2,n}^{2}$ where

Take the union $\widehat{I}_{n}:=\widehat{I}_{1,n}\bigcup\widehat{I}_{2,n}$. Perform spectral decomposition of $\mathbf{S}_{\widehat{I}_{n},\widehat{I}_{n}}$ . Estimate $\theta_{\nu}$ by augmenting the $\nu$-th eigenvector, with zeros in the coordinates $\{1,\ldots,N\}\setminus\widehat{I}_{n}$, for $\nu=1,\ldots,\widehat{M}$. Call this vector $\widehat{\theta}_{\nu}$.

Perform a coordinatewise “hard” thresholding of $\widehat{\theta}_{\nu}$ at threshold

and then normalize the thresholded vectors to get the final estimate $\overline{\theta}_{\nu}$.

Remark : The scheme is specified except for the “tuning parameters” $\gamma_{1}$,$\gamma_{2}$,$\gamma_{3}$ and $\kappa$. The choice of $\gamma_{i}$’s is discussed in the context of deriving upper bounds on the risk of the estimator. It will be shown that, it suffices to take $\gamma_{1}=4$, $\kappa=2+\epsilon$ for a small $\epsilon>0$, and $\gamma_{2}=\sqrt{\frac{3}{2}}\kappa$. An analysis of the thresholding scheme is not done here, but in practice $\gamma_{3}=3$ works well enough, and some calculations suggest that $\gamma_{3}=2$ suffices asymptotically.

3 Estimation of MM

Let $\overline{\gamma}_{1},\gamma_{1}^{\prime}>0$ be such that $\overline{\gamma}_{1}>\gamma_{1}^{\prime}$. Define

The choice of $\gamma_{1}^{\prime}$ and $\overline{\gamma}_{1}$ is discussed in Section 8.5.

Remark : Sparsity of the eigenvectors is an implicit assumption for ASPCA scheme. However, in practice, and specifically with only moderately large samples, it is not always the case that ASPCA is able to select the significant coordinates. More importantly, the scheme produces a bona fide estimator only when $\widehat{I}_{1,n}$ is non-empty. If this is not the case, then one may use the $\nu$-th eigenvector of $\mathbf{S}$ as the estimator of $\theta_{\nu}$. However, determination of $M$ in this situation is a difficult issue, and without recourse to additional information, one may set $\widehat{M}=0$.

Rates of convergence

In this section we describe the asymptotic risk of ASPCA estimators under some regularity conditions. The risk is analyzed under the loss function (3), and it is assumed that condition BA of Section 2 holds. Further, the parameter space for $\theta=[\theta_{1}:\ldots:\theta_{M}]$, over which the risk is maximized, is taken to be $\Theta_{q}^{M}(C_{1},\ldots,C_{M})$ defined through (10) in Section 3.2, where $0<q<2$ and $C_{1},\ldots,C_{M}>1$.

The following conditions are imposed on the “hyperparameters” of the parameter space $\Theta_{q}(C_{1},\ldots,C_{M})$. Suppose that $\rho_{1},\ldots,\rho_{M}$ are as in C1 given below. Define

Observe that, since $C_{\nu}\geq 1$ for all $\nu=1,\ldots,M$, $\rho_{q}(C)\geq\sum_{\nu=1}^{M}\rho_{\nu}^{q/2}\geq 1$.

$\lambda_{1},\ldots,\lambda_{M}$ are such that, as $n\to\infty$, $\frac{\lambda_{\nu}}{\lambda_{1}}\to\rho_{\nu}$ where $1\equiv\rho_{1}>\rho_{2}>\ldots>\rho_{M}$.

$\log N\asymp\log n$ and $\frac{(\log n)^{2}}{n\lambda_{1}^{2}}\to 0$ as $n\to\infty$.

$\frac{\rho_{q}(C)(\log N)^{1/2-q/4}}{\lambda_{1}^{1-q/2}n^{1/2-q/4}}\to 0$ as $n\to\infty$.

We discuss briefly the importance of these conditions. C1 is a repetition of L1. C2 is a convenient and very mild technical assumption that should hold in most practical situations. Second part of ${\bf C2}$ is non-trivial only when $\lambda_{1}\to 0$ as $n\to\infty$. C3 requires some explanation. It will become increasingly clear that, in order to get a uniformly consistent estimate of the eigenvectors from the preliminary SPCA step, one needs C3 to hold. Indeed, the sequence described in C3 has the same asymptotic order as a common upper bound for the rate of convergence of the supremum risk of the SPCA estimators of all the $\theta_{\nu}$’s. So, the implication is that if C3 holds then the SPCA scheme of Johnstone and Lu (2004) gives consistent estimates.

Remark : Note that, $\frac{1}{nh(\lambda)}\leq\frac{1+c}{n\lambda^{2}}$ if $\lambda\in(0,c)$ and $\frac{1}{nh(\lambda)}\leq\frac{1}{\eta(c)n\lambda}$ if $\lambda\geq c$, for any $c>0$. Since $\rho_{q}(C)\geq 1$, C3 guarantees that

In fact, if $\liminf_{n\to\infty}\lambda_{1}\geq c>0$, then the upper bound in (23) can be replaced by $o((\frac{\log N}{n})^{1/2-q/4})$. It will be shown that this is a common (and near-optimal) upper bound on the rate of convergence of the ASPCA estimate of $\theta_{\nu}$’s. If one compares this with the lower bound given by Theorem 2, it is conjectured that (23) should also be a sufficient condition for establishing that the lower bound defined through (15) is also the upper bound on the minimax risk, at the level of rates. However, since our method depends on finding a preliminary consistent estimator of the eigenvectors (in our case SPCA), the somewhat stronger condition C3 becomes necessary to establish rates of convergence of the ASPCA estimator.

2 Statement of the result

Now we state the main result of this section. The asymptotic analysis of risk is conducted only for the estimator $\widehat{\theta}_{\nu}$ for eigenvector $\theta_{\nu}$, and not for the thresholding estimator $\widetilde{\theta}_{\nu}$. Derivation of the results for $\widetilde{\theta}_{\nu}$ requires additional technical work, but can be carried out. It can be shown that in certain circumstances the latter has a slightly better asymptotic risk property. In practice, the thresholding estimator seems to work better when the eigenvalues are well-separated. The following theorem describes the asymptotic behavior of the risk of the ASPCA estimator $\widehat{\theta}_{\nu}$ under the loss function $L$ defined through (3). $g(\cdot,\cdot)$ is defined by (11).

Theorem 3: Assume that BA and conditions C1-C3 hold. Then, there are constants $K:=K(q,\gamma_{1},\gamma_{2},\kappa)$ and $K^{\prime}:=K^{\prime}(q,M,\gamma_{1},\gamma_{2},\kappa)$ such that, as $n\to\infty$, for all $\nu=1,\ldots,M$,

Remark : The expression in the upper bound is somewhat cumbersome, but the significance of each of the terms in (24) will become clear in the course of the proof. However, notice that, if the parameters $C_{1},\ldots,C_{M}$ of the space $\Theta_{q}(M)(C_{1},\ldots,C_{M})$ are such that,

then, Theorem 3 and Theorem 2 together imply that, under conditions BA, C1-C3, A1 and the condition on the hyperparameters given by (15), the ASPCA estimator $\widehat{\theta}_{\nu}$ has the optimal rate of convergence. The condition (25) is satisfied in particular if $C_{1},\ldots,C_{M}$ are all bounded above.

Remark : It is instructive to compare the asymptotic supremum risk of ASPCA with that of OPCA (or usual PCA based) estimator of $\theta_{\nu}$. A closer inspection of the proof reveals that, if for all sufficiently large $n$,

then for some constant $K^{\prime\prime}>0$, under BA and C1-C3, one can replace the upper bound in (24) by

for some constant $\overline{K}$. This rate is greater than that of OPCA estimator by a factor of at most $\log(n\vee N)$. However, observe that, the bound on the risk of OPCA estimator holds under weaker conditions. In particular, Theorem 1 does not assume any particular structure for the eigenvectors.

Proof of Theorem 2

The proof requires a closer look at the geometry of the parameter space, in order to obtain good finite dimensional subproblems that can then be used as inputs to the general machinery, to come up with the final expressions.

A key tool for our proof the lower bound on the minimax risk is Fano’s lemma. Thus, it is necessary to derive a general expression for the Kullback-Leibler discrepancy between the probability distributions described by two separate parameter values.

Proposition 1: Let $\theta^{(j)}=[\theta_{1}^{(j)}:\ldots:\theta_{M}^{(j)}]$, $j=1,2$ be two parameters. Let $\Sigma_{(j)}$ denote the matrix given by (1) with $\theta=\theta^{(j)}$ (and $\sigma=1$). Let $P_{j}$ denote the joint probability distribution of $n$ i.i.d. observations from $N(0,\Sigma_{(j)})$. Then the Kullback-Leibler discrepancy of $P_{2}$ from $P_{1}$, to be denoted by $K_{1,2}:=K(\theta^{(1)},\theta^{(2)})$, is given by

2 Use of Fano’s lemma

We outline the general approach pursued in the rest of this section. The idea is to bound the supremum of the risk on the entire parameter space by the maximum risk over a finite subset of it, and then to use some variant of Fano’s lemma to provide a lower bound for the latter quantity.

Thus, the goal is to find an appropriate finite subset ${\cal F}_{0}$ of $\Theta_{q}^{M}(C_{1},\ldots,C_{M})$, such that the following properties hold.

If $\theta^{(1)},\theta^{(2)}\in{\cal F}_{0}$, then $L(\theta_{\nu}^{(1)},\theta_{\nu}^{(2)})\geq 4\delta$, for some $\delta>0$ (to be chosen). This property will be referred to as “$4\delta$-distinguishability in $\theta_{\nu}$”.

The element $\theta\in{\cal F}_{0}$ is a unique representative of the equivalence class $[\theta]$, where $[\theta]$ is defined to be the class of $N\times M$ matrices whose $\nu$-th column is either $\theta_{\nu}$ or $-\theta_{\nu}$.

Subject to (1), the quantity $\sup_{i\neq j:~\theta^{(i)},\theta^{(j)}\in{\cal F}_{0}}K(\theta^{(i)},\theta^{(j)})+K(\theta^{(j)},\theta^{(i)})$ is as small as possible.

Given any estimator $\widehat{\theta}$ of $\theta$, based on data $\mathbf{X}_{n}=(X_{1},\ldots,X_{n})$, define a new estimator $\phi(\mathbf{X}_{n})$ (an $N\times M$ matrix) as $\phi(\mathbf{X}_{n})=\theta^{*}$ if $\theta^{*}=\arg\min_{\theta\in{\cal F}_{0}}L(\theta_{\nu},\widehat{\theta}_{\nu})$, where $\widehat{\theta}_{\nu}$ is the $\nu$-th column of $\widehat{\theta}$ (i.e., estimate of $\theta_{\nu}$). Then, by Chebyshev’s inequality,

The last inequality is because, if $L(\theta_{\nu}^{(j)},\widehat{\theta}_{\nu})<\delta$ for any $\theta^{(j)}\in{\cal F}_{0}$, then by the “$4\delta$-distinguishability in $\theta_{\nu}$” (property (1) above), it follows that $[\phi_{\nu}(\mathbf{X}_{n})]=[\theta_{\nu}^{(j)}]$, and hence $[\phi(\mathbf{X}_{n})]=[\theta^{(j)}]$.

Two versions of Fano’s lemma are found to be useful in this context. The following version, due to Birgé (2001), of a result of Yang and Barron (1999) (p.1570-71), is most suitable when ${\cal F}_{0}$ can be chosen to be large.

Lemma 1: Let $\{P_{\theta}:\theta\in\Theta\}$ be a family of probability distributions on a common measurable space, where $\Theta$ is an arbitrary parameter space. Suppose that a loss function for the estimation problem is given by $L^{\prime}(\theta,\theta^{\prime})=\mathbf{1}_{\theta\neq\theta^{\prime}}$. Define the minimax risk over $\Theta$ by

where $T$ denotes an arbitrary estimator of $\theta$ with values in $\Theta$. Then for any finite subset ${\cal F}$ of $\Theta$, with elements $\theta_{1},\ldots,\theta_{J}$ where $J=|{\cal F}|$,

To use Lemma 1 choose $P_{i}$ to be $P_{\Sigma_{(i)}}\equiv P_{\theta^{(i)}}:=N^{\otimes n}(0,\Sigma_{(i)})$, where $\Sigma_{(i)}$ is the matrix $\sum_{\nu=1}^{M}\lambda_{\nu}\theta_{\nu}^{(i)}{\theta_{\nu}^{(i)}}^{T}+I$, and $\theta^{(i)}\in{\cal F}_{0}$ $i=1,\ldots,|{\cal F}_{0}|$, are the distinct values of parameter $\theta$ that constitute the set ${\cal F}_{0}$. Then set $Q_{0}=P_{\theta^{(0)}}$, for some appropriately chosen $\theta^{(0)}\in\Theta_{q}^{M}(C_{1},\ldots,C_{M})$ such that the following condition is satisfied.

where the notation “$\asymp$” means that the both sides are are within constant multiples of each other. Then it follows from (28) and Lemma 1 that,

To complete the picture it is desirable that

A different version of Fano’s lemma, due to Birgé (2001), is needed when ${\cal F}_{0}$ consists of only two elements $\theta^{(1)}$ and $\theta^{(2)}$, so that the classification problem reduces to a test of hypothesis of $P_{1}$ against $P_{2}$.

Lemma 2: Let $\alpha_{T}$ and $\beta_{T}$ denote respectively the Type I and Type II errors associated with an arbitrary test $T$ between the two simple hypotheses $P_{1}$ and $P_{2}$. Define, $\pi_{mis}=\inf_{T}(\alpha_{T}+\beta_{T})$, where the infimum is taken over all test procedures.

3 Geometry of the parameter space

We view the space $\Theta_{q}(C)$, for $0<q<2$, as the $N$-dimensional unit sphere centered at the origin, from which some parts have been chopped off, symmetrically in each coordinate, such that there is some portion left at each pole (i.e., a point of the form $(0,\ldots,0,\pm 1,0,\ldots,0)$, where the non-zero term appears only once). In this connection, we define an object that is central to the proof of Theorem 3.3.

Of course, if $C^{q}\geq m^{1-q/2}$ (or $m_{C}\geq m$ ) then as a convention, $r_{m}(C)=1$. Condition (35) ensures that all the points lying on an $(N,m,r)$ polar sphere such that $r\in(0,r_{m}(C))$, are inside $\Theta_{q}(C)$.

4 A common recipe for Part (a) and Part (b)

In the proof of Part (a) and Part (b) of the theorem, there is a common theme in the construction of ${\cal F}_{0}$. Let $\mathbf{e}_{\mu}$ denote the $N$-vector whose $\mu$-th coordinate is 1 and rest are all zero. In either case, if $\{\theta^{(j)},j=1,\ldots,|{\cal F}_{0}|\}$ is an enumeration of the elements of ${\cal F}_{0}$, then the following are true.

There is an $N\times M$ matrix $\theta^{(0)}$, such that $\theta_{\nu}^{(0)}=\mathbf{e}_{\nu}$.

$\theta_{\mu}^{(j)}=\mathbf{e}_{\mu}$ for $\mu=1,\ldots,\nu-1,\nu+1,\ldots,M$, for all $j=0,1,\ldots,|{\cal F}_{0}|$.

$\theta_{\nu}^{(j)}\in{\cal S}(N,m,r,\nu,J)$ for some $m$, $r$ and $J$. $m$ and $r$ are fixed for all $1\leq j\leq|{\cal F}_{0}|$, but $J$ may be different for different $j$, depending on the situation.

The $\theta^{(0)}$ in (F1) is the same $\theta^{(0)}$ appearing in (31). Also, (26) simplifies to

Moreover, in either case, the points $\theta^{(j)}$ are so chosen that

In other words, the set ${\cal F}_{0}$ is $r^{2}$ distinguishable in $\theta_{\nu}$.

5 Proof of Part (a)

Construct ${\cal F}_{0}$ satisfying (F1)-(F3), with

where $r\in(0,1)$ is such that $(1-r^{2})^{q/2}+r^{q}\leq C_{\nu}^{q}$. Thus, $|{\cal F}_{0}|=N-M$. Verify that (37) holds, in fact the lower bound is $2r^{2}$, with an equality. Therefore, (31) applies, with $\delta=\frac{r^{2}}{2}$. Since $nh(\lambda_{\nu})$ is bounded above, and $\log(N-M)\to\infty$ as $n\to\infty$, (12) follows from (36).

6 Connection to “Sphere packing”

Our proof of Part (b) of Theorem 2 depends crucially on the following construction due to Zong (1999).

7 Proof of Part (b)

Structures of ${\cal F}_{0}$ for the three cases in (14) are similar. Set $m\leq(N-M)$, large. Set $c_{1}=\log(9/8)$, $A_{q}=(9c_{1}/2)^{1-q/2}$. Choose $r\approx\sqrt{\delta_{n}}$, and define the set ${\cal F}_{0}$ satisfying (F1)-(F3) and the following construction.

Set $|{\cal F}_{0}|=|Y_{m}^{*}|$, where $Y_{m}^{*}$ is the set defined in Section 6.6. Set,

where $\mathbf{z}^{(j)}=(z_{1}^{(j)},\ldots,z_{m}^{(j)})$, $j\geq 1$, is an enumeration of the elements of $Y_{m}^{*}$. Observe that, for all $j\geq 1$,

where $supp(\mathbf{z}^{(j)})$ is the set of nonzero coordinates of $\mathbf{z}^{(j)}$. Therefore, (37) and (36) hold for all $j\geq 1$.

Take $m=[nh(\lambda_{\nu})]$ and $r^{2}=c_{1}$. Observe that, for all $j\geq 1$,

Thus, ${\cal F}_{0}\subset\Theta_{q}^{M}(C_{1},\ldots,C_{M})$. Further, since $nh(\lambda_{\nu})\to\infty$, $\log|{\cal F}_{0}|\geq c_{1}nh(\lambda_{\nu})(1+o(1))$. Since (37) and (36) hold, with $\delta=\frac{r^{2}}{4}$, from (31) the result follows, because

Take $m=N-M$ and $r^{2}=\frac{c_{1}(N-M)}{nh(\lambda_{\nu})}$. Then, for all $j\geq 1$,

The result follows by arguments similar to those used for the case $nh(\lambda_{\nu})\leq\min\{c_{1}(N-M),A_{q}\overline{C}_{\nu}^{q}(nh(\lambda_{\nu}))^{q/2}\}$.

Take $m=[c_{1}^{-q/2}(9/2)^{1-q/2}\overline{C}_{\nu}^{q}(nh(\lambda_{\nu}))^{q/2}]$ and $r^{2}=c_{1}\frac{m}{nh(\lambda_{\nu})}$. Again, verify that $m\to\infty$ as $n\to\infty$ (by A1), and for $j\geq 1$,

and the result follows by familiar arguments.

7.4 Proof of (15)

The construction in all three previous cases assumes that the set of non-zero coordinates is held fixed (in our case $\{M+1,\ldots,M+m\}$) for every fixed $m$. However, it is possible to get a bigger set ${\cal F}_{0}$ satisfying the requirements, if this condition is relaxed.

Suppose that $A_{q,\alpha}=(\alpha/2)^{1-q/2}$, and the condition in (15) holds for some $\alpha\in(0,1)$. Set $m=[(\alpha/9)^{-q/2}(9/2)^{1-q/2}\overline{C}_{\nu}^{q}(nh(\lambda_{\nu}))^{q/2}(\log N)^{-q/2}]$ and $r^{2}=(\alpha/9)\frac{m}{nh(\lambda_{\nu})}$. Take $c_{q}(\alpha)=(\alpha/9)^{1-q/2}$. Observe that $m\to\infty$ as $n\to\infty$, $m=O(N^{1-\alpha})$ and $r\in(0,1)$. Set $\theta^{(0)}=[\mathbf{e}_{1}:\ldots:\mathbf{e}_{M}]$. For every set $\pi\subset\{M+1,\ldots,N\}$ of size $m$, construct ${\cal F}_{\pi}$ satisfying (F1)-(F3) such that,

As before, ${\cal F}_{\pi}\subset\Theta_{q}^{M}(C_{1},\ldots,C_{M})$, for all $\pi$, so that (36) and (37) are satisfied. Let ${\cal P}$ to be a collection of such sets $\pi$ such that, for any two sets $\pi$ and $\pi^{\prime}$ in ${\cal P}$, the set $\pi\cap\pi^{\prime}$ has cardinality at most $\frac{m_{0}}{2}$. This ensures that

This also ensures that the sets ${\cal F}_{\pi}$ are disjoint for $\pi\neq\pi^{\prime}$, since each $\theta_{\nu}^{(j)}$ for $\theta^{(j)}\in{\cal F}_{0}$ is nonzero in exactly $m_{0}+1$ coordinates. Define ${\cal F}_{0}=\bigcup_{\pi\in\cal P}{\cal F}_{\pi}$. Then

By Lemma 7, stated in Section 9.4, there is a collection ${\cal P}$ such that $|{\cal P}|$ is at least $\exp([N{\cal E}(m/9N)-2m{\cal E}(1/9)](1+o(1)))$, where ${\cal E}(x)$ is the Shannon entropy function :

Since ${\cal E}(x)\sim-x\log x$ when $x\to 0+$, it follows from (42) that,

since $m=O(N^{1-\alpha})$. Finally, observe that

8 Proof of Part (c)

Consider first the proof of (16). Fix a $\mu\in\{1,\ldots,M\}\setminus\{\nu\}$. Define $\theta^{(1)}$ and $\theta^{(2)}$ as follows. Set $r^{2}=\frac{2}{ng(\lambda_{1},\lambda_{2})}$ (assume w.l.o.g. that $r<1\wedge C_{0}$). Take $\theta_{\mu^{\prime}}^{(j)}=\mathbf{e}_{\mu^{\prime}}$, $j=1,2$ for all $\mu^{\prime}\neq\mu,\nu$. Define

Observe that $\theta_{\mu}^{(j)}\perp\theta_{\nu}^{(j)}$, $j=1,2$, $\langle\theta_{\nu}^{(1)},\theta_{\nu}^{(2)}\rangle=\sqrt{1-r^{2}}=\langle\theta_{\mu}^{(1)},\theta_{\mu}^{(2)}\rangle$ and $\langle\theta_{\mu}^{(1)},\theta_{\nu}^{(2)}\rangle=r=-\langle\theta_{\nu}^{(1)},\theta_{\mu}^{(2)}\rangle$. Also, by A1, $\theta^{(j)}\in\Theta_{q}^{M}(C_{1},\ldots,C_{M})$, for $j=1,2$.

Let $P_{j}=N^{\otimes n}(0,\Sigma_{(j)})$. Then

Apply Lemma 2 for testing $P_{1}$ against $P_{2}$. Define $p_{mis}=\inf_{T}(\alpha_{T}\vee\beta_{T})$ and observe that $p_{mis}\leq\pi_{mis}\leq 2p_{mis}$. Since the lower bound in (33) is symmetric w.r.t. $\pi_{mis}$, and $\pi_{mis}$ is symmetric w.r.t. $P_{1}$ and $P_{2}$, it follows that

Since, $L(\theta^{(1)},\theta^{(2)})=2(1-\sqrt{1-r^{2}})\geq r^{2}$, and $r^{2}=\frac{2}{ng(\lambda_{\mu},\lambda_{\nu})}$, use (28) with ${\cal F}_{0}=\{\theta^{(1)},\theta^{(2)}\}$ and $\delta=r^{2}$ to get,

Now, let $\mu$ vary over all the indices $1,\ldots,\nu-1,\nu+1,\ldots,M$ and the result follows.

In the situation where $\overline{\delta}_{n}\not\to 0$, as $n\to\infty$, simply take $\mu$ ($\neq\nu$) to be the index for which $g(\lambda_{\mu},\lambda_{\nu})$ is minimum. Then apply the same procedure as in above with $r\in(0,C_{0})$ fixed.

Proof of Theorem 1

We require two main tools in the proof of Theorem 1 - one (Lemma 5) is concerned with the deviations of the extreme eigenvalues of a Wishart$(N,n)$ matrix and the other (Lemma 6) relates to the change in the eigen-structure of a symmetric matrix caused by a small, additive perturbation. Sections 9.1 and 9.2 are devoted to them. The importance of Lemma 6 is that, in order to bound the risk of an estimator of $\theta_{\nu}$ one only needs to compute the expectation of squared norm of a quantity that is linear in $\mathbf{S}$ (or a submatrix of this, in case of ASPCA estimator). The second bound in (132) then ensures that the remainder is necessarily of smaller order of magnitude. This fact is used explicitly in deriving (66).

Remark : In view of Lemma 6, $H_{\nu}(\Sigma)$ becomes a key quantity in the analysis of the risk of any estimator of $\theta_{\nu}$. Observe that,