Classification of high-dimensional data with spiked covariance matrix structure
Authors: Yin-Jen Chen
Minh Tang
Venue: Accepted by TMLR ()

Abstract: We study the classification problem for
high-dimensional data with $n$ observations on $p$ features where the
$p \times p$ covariance matrix $\Sigma$ exhibits a spiked eigenvalue structure and the
vector $\zeta$, given by the difference between the {\em whitened} mean
vectors, is sparse. We analyze an adaptive
classifier (adaptive with respect to the sparsity $s$) that first
performs dimension reduction on the feature vectors prior to classification in
the dimensionally reduced space, i.e., the classifier whitens
the data, then screens the features by keeping only those corresponding
to the $s$ largest coordinates of $\zeta$ and finally applies Fisher
linear discriminant on the selected features. Leveraging recent
results on entrywise matrix perturbation bounds for covariance
matrices, we show that the resulting classifier is Bayes optimal
whenever $n \rightarrow \infty$ and $s \sqrt{n^{-1} \ln p}  \rightarrow
0$. Notably, our theory also guarantees Bayes optimality for the corresponding quadratic discriminant analysis (QDA). Experimental results on real and synthetic data further indicate that the proposed approach is competitive with state-of-the-art methods while operating on a substantially lower-dimensional representation.

Published in Transactions on Machine Learning Research (January 22026)
Classification of high-dimensional data with spiked covari- ance matrix structure
Yin-Jen Chen
acgooguorkk@maail.com com
Deportment of Statistics North Carolina State University
Minh Tang
miangg8@nccuu edu
Department Statistics
North Carolina State University
Reviewed on OpenReview; https: opw.rvvieee.......... net/ fttt/f/rmm id= ToQDtTaaas
Abstract
We study the classification problem for high-dimensional data with n observations on features where the p x p covariance matrix E exhibits a spiked cigenvalue structure and the vector (; given by the difference between the whitened mean vectors, is sparse. We analyze an adaptive classifier (adaptive with respect to the sparsity s) that first performs dimension reduction the feature vectors prior to classification in the dimensionally reduced space, ie., the classifier whitens the data, then screens the features by keeping only those corre- sponding to the largest coordinates of and finally applics Fisher linear discriminant on the selected features. Leveraging recent results on entrywise matrix perturbation bounds for covariance matrices, we show that the resulting classifier is Bayes optimal whenever - oo and In O. Notably, our theory also guarantees Bayes optimality for the corresponding quadratic discriminant analysis (QDA). Experimental results on real and syn- thetic data further indicate that the proposed approach is competitive with state-of-the-arl methods while operating on substantially dinee-dimesioonll representation.
1 Introduction
Classification is one of the most important and widely studied inference tasks in statistics and machine learn- ing. Among standard classifiers, the Fisher lincar discriminant analysis (LDA) rule is ospecially popular for its ease of implementation and interpretation. More specifically, suppose that we are given a p-variate ran- dom vector drawn from a mixture of two multivariate normal distributions NINa(Hi (1- 2) and our goal classify Z into one of the two classes. The Fisher LDA rule is then given by
1 if (I- -)))a-m))]]000) (Z #1 In
(z)
I-w
(1.1)
2 if (E- (M2- m1)) (Z 10 In
provided that Σ positive definite. The Fisher rule TF the Bayes decision rule, i.c., it achieves the smallest. mis-classification error with respect to 0-1 loss for classifying Z. Tr is, however, not directly applicable in practice it involves the unknown parameters Σ-1 ln and μ2; we thus usually compute the sample covariance matrix Σ and the sample means X, and X, from given training data set with n observations and then plugged these quantities into Eq. I. 0 the resulting classifier is termed as the plug-in LDA rule r-
The classification accuracy of i well-understood the low-dimensional regime where p & n. In particular fr is asymptotically optimal, i.c., achieves the Bayes error rate as 3o for fixed p. This behavior, however, might longer holds in the high-dimensional settings p E p where the sample

---

Published in Transactions on Machine Learning Research (anuary/20200
covariance matrix Σ is singular. Indeed, Bickel & Levina (2004) proved that the error rate for TF where we replace Σ-1 with the Hoore-Penross, pseudo-inverse of X could be as bad as random guessing and that, furthermore, the naive Bayes (NB) rule which ignores the correlation structure in > typically outperforms
F...
Continuing this line of inquiry, other independence rules (TR) have been proposed in Tibshirani et al. (2002) and Fan & Fan (2008) but these classifiers have two potentially major drawbacks, namely (I) the accumula- tion of error in estimating μi and A22 and (2) mis-specification of the covariance matrix Σ. In particular Fan & Fan (2008) demonstrated that the noise accumulation in estimating μi and μ2 alone is sufficient to degrade the performance of IR classifiers, Shao et al. (2011) extended this result to show that the noisc accumulation can also lead to the plug-in LDA classifier being, asymptotically, no better than random guessing even when E is known, and finally Fan et al.) (2012) showed that ignoring the correlations structure in Σ prevents IR. classifiers from achieving Bayes optimality. These results indicate the need for imposing sparsity conditions on the features and the important role of feature selection for mitigating the estimation crrors associated with growing dimensionality.
To address the above limitations of IR classifiers, Shao et al. 2011 proposed thresholding of both Σ and
#22 M11 this is similar to the motivation for regularized covariance estimators in Bickel & Levina (2009)
In contrast, Cai & Liu (2011); Fan et a (2012) 2): Mai et al.] (2012): : Witten & Tibshirani (2011) and Cai &
Zhang 2019 imposed sparsity conditions on the discriminant direction B == Σ- 'p- - -11 and use penalized
estimation approaches to recover B. Caa & Liu (2011) also noted two potential advantages to this approach, namely that the assumption of sparsity on 3 less restrictive than assuming sparsity for both Σ (or E and /12 and secondly Yr only depends on E-1 and /22 141 through their product (12- and thus consistent estimation of 3 is sufficient. The sparsity assumption on 8 leads to procedures and results that resemble those for high-dimensions formulated in terms of a linear model. For example, Cai & Liu (2011) and Cai & Zhang 2019)) considered a lincar programming approach similar to the Dantzig selector (Candes & Tao 2005 while Mai et al. (2012) studied a sparse discriminant analysis rule that used the Lasso (Tibshirani] 1996). Nevertheless it had been observed that, empirically, these approaches can lead to classification rules which select a larger number of features than necessary, and one possible explanation is that the correlation structure in E also induced correlations among the entries for any estimate / of B.
linear regression even though the classification problem is generally not
In this paper we consider a different approach where we first perform dimension reduction on the feature vectors (using PCA) prior to classification (using LDA) in the dimensionally reduced space. As PCA is an important and ubiquitous pre-processing step in high- sigh-dimensional data analysis, there is a sizable number of work devoted to this approach. We refer to a generic classifier from this combination as Ida o pea. For example Section 9.1 of Jolliffe (2002) provides detailed review of combining LDA with different variants of PCA in the low-dimensional setting while NNu et al. (2015) proposed the use of the reduced rank LDA together with class-conditional PCA in the high- sigh-dimensional settings. In terms of applications, Ida o pca is also used for faces and images recognition (Zhao et al. 1998 Belhumeur et al. 1997 Prasad et al. 2010
and recovering genetic patterns (Jombart et al. 2010)
While Ida o pca classifiers arise quite naturally, their theoretical properties in the high-dimensional setting remains an open problem. In particular their analysis requires possibly different techniques and assumptions compared to those based on direct estimation of the discriminant direction 33. More specifically consistency results for 3 are usually based ideas from high-dimensional regression including assumptions on bounded and concentrated eigenvalues of Σ, see e.g., Bickel & Levina 22004); Fan & Fan (2008); Shao et al. (2011)
Cai & Zhang (2019). However Wang & Fan (2017) noted that bounded eigenvalues are incompatible with the presence of strong signals (eigenvalues) in the data and might. be problematic in fields such as genomics,
economics and finance. These assumptions are nevertheless imposed mainly due to the limitation of quanti- fying the estimation error for Σ in terms of the spectral norm difference for -ー Σ. In contrast the idealized setting for dimension reduction via PCA is when Σ contains a small subset of signal eigenvalues that accounts for most of the variability in Σ. This idealized setting for PCA is also distinct from the idcalized setting for graphical models wherein E-1 is typically assumed to be sparse.
2

---

Published in Transactions on Machine Learning Research (January /202)
In summary, our contributions in this paper follows. We analyze the theoretical properties of prototypical Ida o pea classifier under a spiked-covariance structure assumption a widely-adopted covari- ance model for high-dimensional data - where E contains a few large eigenvalues that are well-separated from the remaining (small) cigenvalues. In particular we show in Section that Ida o pea is asymptotically Bayes-oplimal as n + 00 and n 1 Inp + This is, to the best of our knowledge, the first Bayes optimal consistency result for classification after performing dimension reduction via. PCA. In Section f we demon- strate empirically, for both simulated and real data, that Ida o pca selects fewer features while also having error rates that are competitive with existing classifiers based on estimating the discriminant direction B. The theoretical and numerical results provide a clear example of the synergy linking dimension reduction with classification. Finally, in Section 5 we extend Ida o pea to classify data with (1) K 2 3 classes or (2) two classes but with unequal covariance matrices or (3) feature vectors that are elliptical dist but not necessarily multivariate normal.
2 Methodology
2.1
Notation and settings
For a vector x E RP, the conventional to quasi-norm and the f &y and Es norms are denoted by ||aloo |<<|] ||~|| and ||olloo respectively. For p N, we denote the set [1.-.. by [p]. Given x R2 and a non-empty set A [p], we write JA (zz: A) to denote the column vector obtained by keeping only the elements in whose indices belong to A. The operation "o", when applied to matrices, represents the Hadamard (entrywise) product. For e [p]. e()) the ith standard basis vector of R2 and 1, is a vector in TR" whose elements are all I; also write e, and when the choice of p is clear from context and Ip stands for the identity matrix in Rexp, For I, E IP, the standard Euclidean inner product between and y is denoted as (zz y) := y z. For a matrix M E RPxq, the Frobenius norm and spectral norm of M are written as || p and || respectively. We will omit the subscript in ||-1|2 when it is clear from context that the spectral norm or l2 norm is intended.
The two-to-infinity norm of M defined as
|M = sup |=1==1
|
max || |M]]
(2.1)
/E|p]
where M; represent the ith row of M. We note that || < MIM < || Let tr(M) and |M denote the trace and determinant of a square matrix M. The effective rank of a square matrix M is defined as r(M) == tr(M)/ |MI]; the effective rank of a matrix is a useful surrogate measure for its complexity, see e.g.: Vershynin 2018). Let O()), o(), O(-) and 2(-) represent the standard big-O, little-o, big-Thete and big- Sig-Omega relationships. Finally, 1 (AA and |A] stand for the indicator function of set A and its cardinality, respectively.
In the subsequent discussion we shall generally assume, unless specified otherwise, that we have access to training set (X.....,, X1n,} and (X2i...., X2ng} whose elements are independently and identically dis- tributed random vector from the p-variate distributions Np Σ) (class 1) and Nn 2) (class 2). The sample means for each class and the pooled sample covariance matrix are denoted as
n2
Xi
X1j+
X₂ =
X2j1
(2.2)
112
j=1
=
X1 - XX)
+ Σ
(X2j Xx) - x2)")
(2.3)
j=1
On numerous occasions we also need a variant of the sample covariance matrix where we replace the sample means with the true means. We denote this matrix as
Eo 1 (Xn, -1)( -T)T (X2s (2) -22) n
n2
(2.4)
3

---

Published in Transactions on Machine Learning Research (January /2026)
2.2 Linear discriminant analysis and whitening matrix
Let Σ be a p x p positive definite matrix. The whitening matrix W is a linear transformation satisfying WWT = E-', and is generally used to decorrelate random variables and scale their variances to 1. The whitening transformatinn unique only up to orthogonal transformation as WTTTWT = WW for any p orthogonal matrix TV see Kessy et al. (2015) for comparison between several common choices of whitening transformation. Hence, for this paper, we take W E the wrique positive semidefinite square root of Σ- Note that --1/2 is the whitening transformation which minimizes the expected mean square error between the original data and the whitened data (Eldar & Oppenheim 2003). Given W, we define the whitened direction as Ç = W(a2 - n1) and let Se = U S + 0}. We referred to the elements of Se as the whitened coordinates or variables. The Fisher linear discriminant rule is then equivalent to
[1 if (w(Z- (())) = C3, (w(z - K Sc
S In
m
Tr(Z)= 2 otherwise,
I-mi
(2.5)
where Z ~ (TNpp(11 Σ) - (1- M1)Np( Σ). Recall that, for vector & and set of indices A, the vector LA is obtained from & by keeping only the elements indexed by A. Since 3 = E-1(22 -11 is the Bayes direction, a significant portion of previous research is devoted to recovering 3 and the discriminative set S3 = {j : By * 0}. We refer to the elements of Sg as the discriminative coordinates or variables. We emphasize that general the discriminative set and whitened set are not the same, except for when E has certain special structures, c.g., E being diagonal.
p < then the empirical whitening matrix W = E-122 is well-defined and the empirical whitened variables Ç W(X, X1) are approximately decorrelated in the plug-in LDA rule, ie., by the law of large numbers, WEW Ip and hence Varla] cI,, for some constant e 0. In contrast, the coordinates of the estimated discriminant direction 3 are not decorrelated since Var[3] xs cx- It is thus easier to quantify the impact of any arbitrary whitened feature toward the classification accuracy than to quantify the impact of an arbitrary raw feature.
2.3 Combining LDA and PCA
Eq. (.5 can be viewed under the framework of first performing dimension reduction and then doing clas- sification in the lower-dinensional spaces with the main focus being that. thresholding of the discriminant direction corresponds to dropping a subset of the raw features while thresholding of the whitened direction corresponds to dropping a subset of the whitened features. We now discuss the estimation of the whitcning matrix W.
If n x p then the sample covariance matrix 3 often- times ill-conditioned, and is furthermore singular when < p. The estimation of either the precision matrix Σ-1 the whitcning matrix W is thus challenging in these regimes, especially since Σ-1 and W both contains O(p?) entries. An universal approach to address this difficulty is to reduce the number of parameters needed for estimating Σ, for example by assuming that E has a parametric form with o(p²) parameters, or by introducing regularization terms to induce certain structures in the estimate E.
The use of regularized estimators for covariance matrices are popular in high-dimensional classification, see e.g. Bickel & Levina (2009); Shao et al.] (2011), and theoretical analysis for these estimators are usually based assumptions about the sparsity of either E or E-1 and then bounding the estimation error in terms of the spectral norm differences || - or | E| There are, however, two potential drawbacks to this approach. Firstly, quantifying the cstimation error in terms of the spectral norm can be quite loose and the resulting bounds might not capture the difference in geometry between the eigenspaces of Σ and that of the perturbed matrix E (e.g.. Johnstone (2001) observed that the first few largest eigenvalues of the sample covariance matrix E are always larger than those for in the high-dimensional sctting). Secondly, while it is not always appropriate assume sparsity of Σ or E-1 -- & Liu (2011) noted that some form of sparsity is needed to guarantee consistent estimation of Σ E-1 under spectral norm error.
---1~>-11
4

---

Published in Transactions on Machine Learning Research (January /202)
In this paper we assume different structure for Σ, namely that E is spiked-covariance matrix with few leading eigenvalues (the spike) that are well-separated from the remaining eigenvalues (the bulk). More specifically, we assume that D satisfies the following condition.
Assumplion 1 (Spiked covariance matrix). Let u1,...; ud be orthonormal vectors in IP and assume that the covariance matrix Σ for the p-variate distributions NG(un, Σ) and (12 Σ) is of the form
M Σ u u + I -UN + I
(2.6)
Here A diag(As) is d' d diagonal matrix, 4= and Zp is the identity matrix. We assumne implicitly that \1 2 ... - Aa 0, > 0 and d <
1=(()) (ux) &s [d] is xx x matrix with orthonormal columns,
Covariance matrices with spiked structures have been studied extensively the high-dimensional statistics literabure, see e-g.: Johnstone (2001) and Chapter 11 Yao et al. (2015), and there is significant number of results for consistent estimation of U under the assumption that the support of U is sparse, e.g., either that the number of non-zero rows of U is small compared to p or that the ly quasi-norm, for some q E [0, 1]. of the columns of U are bounded. The case when g = 0 and g > 0 correspond to "hard" and "soft" " sparsity constraints, respectively; see for example Birnbaum et al. (2012); Berthet & Rigollet (2012); Vu & Lei (2012ps Cai et al. (2013) and the references therein. In contrasts to the above cited results, in this paper we do not impose sparsity conditions on U but instead assume that U have bounded coherence, namely that the maximum lg norm of the rows of U are of order O(p-1/2, see Assumption for a precise statement. The resulting matrix Σ will no longer be sparse. The main rationale for assuming bounded coherence is that the spiked eigenvalues A can grow linearly with p while still guaranteeing bounded variance in Σ. There is thus large gap between the spiked eigenvalues and the bulk eigenvalues, and this justifies the use of PCA as a pre-processing step.
If Σ satisfies Condition ! then the whitening matrix is given by
1 -I
W WDDDu
(2.7)
where D = diag(nk) is d x d diagonal matrix with diagonal entries mh = (Ax + 0 The above form for W suggests the following estimation procedure.
1. Extract the d largest eigenvalues and corresponding cigenvectors of the pooled sample covariance matrix E. Let A denote the diagonal matrix of these d largest eigenvalues and ú denote the x d orthogonal matrix whose columns are the corresponding eigenvectors.
2. Estimate the non-spiked eigenvalues by
0² tr(E) -r((A)
(2.8)
p-6
3. Let = (A &II)) and estimate W by w = ÜDūT + +0- -ûūT).
A Although the sample covariance matrix E is generally a poor estimate of Σ, when n < p, the eigenvectors L7 corresponding to the d largest cigenvalues of E are neverthcless accurate estimates Fan al. (2018) and Cape minimum is taken over all d x d orthogonal transformations T;; also Theorem [ of the current paper. We can thus transform Á by an orthogonal transformation T so that the resulting (transformed) rows of ú are uniformly close to the corresponding rows of 24
U. In particular, al. (2019) provide uniform error bound for min7 - U||2--x where the
Given the above estimate for W, we then have the prototypical Idaopea classifier in Algorithm Algorithm [ contains two tuning parameters, namely d, the number of principal components in the PCA step and s, the number of coordinates of estimated whitened direction & that we preserve. The choices for d and s correspond to the number of spiked eigenvalues in Σ and the sparsity level of the whitened direction Ç. We show in
5

---

Published in Transactions on Machine Learning Research (January 2002)
Section 8 that, under certain mild conditions, one can consistently estimate using Eq. 336 Consistent estimation of d can be obtained using results in Bai & Ng (2002); ) Alessi et al. (2010) Hallin & Liška (2007)
among others. Algorithm [ thus yields a classifier that is adaptive with respect to both d and s.
Algorithm 1: Ida pea decision rule
Input: X1, X2, E and the test sample Z
Output: Tidaopca(Z)
Algorithm
// Step 1: Perform PCA on the feature vectors for the training data (standard PCA approach)
// Step 2: Extract the d largest principal components and obtain W ÜDū + -ûú
// Step 3: Take x; = WX, i = 1,2 and X, = 0.5(X, + X2). Let Ç= X2 - X1 and form the
indices set & by selecting the largest coordinates of Č in modulus.
// Step 4: Given the test sample Z, take Z WZ and plug the sub-vector of ?, and X. corresponding to the indices in & into the Fisher discriminant rule, i.e.
if c 2 I In 22
1 Idaopea (Z) =
(2.9)
otherwise
2.4
Related works
2.4.1
Whitening matrix and spiked covariance
We first discuss the relationship between Algorithm I and two other relevant classifiers, namely the features annealed independence rule (FATR) of & Fan 2008 and the LDA with CAT scores (CAT-LDA) of Zuber & Strimmer (2009). FAIR is the independence rule applied to the features prescreened by the two-sample test while CAT-LDA scores decorrelate the statistios by first whitening the data using the sample correlation matrix diag(2 ) 1/12 & diag(E)-1/2 CAT-LDA is motivated by the empirical observation that accounting for correlations is essential in the analysis of proteomic and metabolic data.
If there are no correlation or if the correlations are negligible, i.e., when is diagonal or approximately diagonal, then Algorithm | is essentially equivalent to both FAIR and CAT-LDA. However, when the cor- relations not negligible, the performance of IR classifiers such as FAIR degrades significantly due to mis-specification of the covariance structure (Fan al. 2012) Both Algorithm [ and CAT-LDA apply whitening transformation before doing LDA but the choice of whitening matrices are different between the two procedures. In particular CAT-LDA uses the asymmetric whitening transformation Corr(E) ///dgg(EE) where Corr(E) is the matrix of correlations. The choice of whitening transformation is not important if Σ is known or if we are in the low-dimensional setting where < n for which consistent estimation of Σ is straightforward. It is only when X n or > n when the different regularization strategies used in CAT-LDA and Algorithm D lead to possibly different behaviors for the resulting classifies; see for example the numerical comparisons given Section f of the current paper. Finally we view the theoretical results in Scction 3 as not only providing justification for Algorithm D but also serve as examples of thcoretical analysis that can be extended to other classifiers in which dimension reduction is done prior to performing LDA. Indeed, to the best of our knowledge, there are no theoretical guarantees for the error rate for the CAT-LDA classifier in the high-dimensional setting.
Finally we note that our paper is not the first work to study higl-dimensional classification under spiked covariance matrices of the form in Condition [ In particular [ifaou et al. (2020a) also analyzed the per- formance of LDA under spiked covariance structures; however, their thcoretical results are premised under different settings and framework from ours and most importantly they do not show that their classification rule achieves the Bayes error rate; Remark i for detailed comparisons of the assumptions and theoretical
6

---

Published in Transactions on Machine Learning Research (annuary/20000
results in our paper against that of Sifaou et al. 2020a With different focus, Hao et al. 2015 examined the effects of spiked eigenvalues when trying to find an orthogonal transformation T of the discriminant
direction / so that T3 is sparse. More specifically they show that if 2 is known then T can be computed using the cigen-decomposition of the matrix Enot where
Eoo = =I = ~
for a given Y 0
(2.10)
where Ap = 122 Pa is the vector of mean differences between the two classes. Let Urot.m be the px matrix whose columns consist of the orthonormal eigenvectors of Eyot corresponding to the m largest eigenvalues. When m = then Mtotp diagonalizes Eno and Hao et al. (2015) showed that Meop further sparsifies / in that 1 Bllo I d+1 whenever satisfies the spiked covariance assumption in Assumption and hence it might be beneficial rotate the data before performing classification.
When E unknown Hao et al. (015)) propose the following procedure for estimating Uno for some choice of m mnn p}: here 7 0 is a user-specified parameter.
Algorithm 2: Rotation preprocessing in LDA
Algorithm
// Step 1: Perform PCA on Êno E (X2 (x) X1)J to extract the m largest principal
components. Let Uroa,m be the resulting p x m matrix.
// Step 2: Rotate the training set to Xii, N....,,,fotnn oot. where m E R for
E [1.2] and j {1,2..... nz).
// Step 3: Apply some discriminant direction based LDA method such as Fan et al. (012)) ; Cai &
Liu 2011
Mai et al. (2015) on the rotated data set (U X
We will use this pre-processing step in some of our numerical experiments Section Nevertheless we emphasize that the thcoretical results of Hao et al. (2015) assume is known as they are mainly concerned with the analysis different approaches for sparsilying B, i.e., they explicitly chose not to address the important issue of how the estimation error for Urorm impacts the classification accuracy.
2.4.2 Connection to high-dimensional sparse LDA
Let X be the (mi + n12) x matrix whose rows are the {Xij) and y E R" be the vector whose first ny elements are set to (n/nn)) and the remaining n2 elements are set to n/n22. Mai et al. (2012) reframed LDA in high-dimension as the solution to a Lasso-type problem (we have omitted the intercept term for simplicity of presentation)
3 = arg min SERP
| X,81 1|
(2.11)
for some À > In a similar spirit Eq. 2.111, we might, conceptually, also reformulate Algorithm as a (general) Lasso-type problem (Tibshirani & Taylor 2011)
i arg min BERP
p X,3 1/11
(2.12)
2
for some À > 0, where V = 1/ (once again omitting the intercept for simplicity). As Σ is invertible,
Eq. (212 is equivalent to first solving
} ID- xc12 + X||lI
= arg min CERP
(2.13)
112 and then setting 8 - 1/2 Comparing Eqs. (2.12) and (2.13) against Eq. 2..1 one could argue that the main difference between lda pea in Algorithm and the lassoed LDA classifier of Mai al. (2012) due to the different choice of the predictor variables X X in the Lasso regression.
where X
7

---

Published in Transactions on Machine Learning Research (January /2026)
This argument, however, overlooks the important fact that the transformation E-1/2 needs to be estimated in Algorithm | while for the general Lassso in Eq. (2.12), the matrix v is specified a priori and thus free from any estimation error. Indeed, while there have been significant efforts devoted to understanding the
solution path of Eq. (.12) in the high-dimensional setting - see, for example, Duan et al. (2016): Tibshirani
& Taylor (2011) p among others- its behavior when is determined empirically from the data is yet to be theoretically investigatod.
Next note that, instead of Eq. 2..1 ,, one can also estimate via an & optimization problem similar to that of the Dantzig selector (Candes & Tao 2005 Cai & Liu 2011 Cai & Zhang 2019 namely
= =rg min [|5ll: subject to | (X2 -1)l s 2]. CERP
(2.14)
The main difference between Algorithm and Eq. (2.14) that Algorithm first whitens the feature vectors using an empirical estimate of the covariance matrix followed by feature selection on the whitened data (see Eq. (2.13)) In contrast, Eq. (2.14) attempts to find Ç for which its transformation Σ1/²5 is most similar to X2 X1. While the approaches underlying Algorithm [ and Eq. (2.14) are quite similar, we believe that Algorithm
1 provides a more direct link between PCA and LDA for high-dimensional classification as it is particularly suitable for high-dimensional data generated from a (possibly low-dimenisonal) factor model (where the covariance matrix will now have a spiked structure with leading eigenvalues that grow with the dimension p). Soe Fan et al. (2021) for further discussion of factor models and their applications to statistics and machine learning.
Although Algorithm ! and Eq. 2.11 share conceptual link in optimization formulation, they differ signif- icantly in algorithmic structure, leading to distinct computational profiles. We now compare the computa- tional complexity of Ida pea and Lassoed LDA. For our method, the dominant cost comes from computing the top d singular values of the centered data matrix X, with cost O(ppd), where d upper-bounds the rank d; see Halko et al. (201I) and Feng & Yu (2023) for more details. In contrast, Lassoed LDA solves an f1-penalized least squares problem with complexity Omn/T), where is the number of iterations until con- vergence. The iteration count T' can be substantial in practice, especially when regularization is weak or X is ill-conditioned as often occurs in settings with low-rank highly correlated features. Consequently, our method provides substantial computational advantages in high-dimenssonall latent variable models, while maintaining competitive classification accuracy and working in a
space (see Section et al. (2024) assumed that the discriminant direction 3 (and the whitening direction Ç as in our work) lies entirely within the column space spanned by the d leading principal components U. However, as noted by Jollifee (1982), the principal components associated with small eigenvalues can be just as important those associated with large eigenvalues in real data analysis (see also the simulation settings for Model 1 and Model 3 in Section 4. By working with the whitened data, Algorithm I also takes into account the non-leading principal components and furthermore, by performing feature selection in the whitened space, avoid the need to specify a priori the low-dimensional subspaces containing B3 and/or Ç.
Finally we conclude this section by comparing Algorithm [ with the principal component classificrs proposed in Bing et al.] (2024); Bing & Wegkamp 2003).. J. In particular both Bing & Wegkamp) (2023) and Bing
3 Theoretical properties
In this section, we derive the theoretical properties of the Ida pea classifier in Algorithm the case where the feature vectors X are sampled from mixture of two multivariate Gaussians. Extensions of these results to the case where X is mixture of K > 2 elliptical, but not necessarily multivariate normal, distributions are discussed in Section 7.3 and Section 5.1
If X ~ (1 Σ) + (1 - 31 N1) Σ) then the Fisher's rule Yr has error rate
R . + |K||| + |<||| In 1 1-11)=))))))) + (1 ) 1 ||||| - IIsll In 1 mi
(3.1)
- Here & is the cumulative distribution function for ,O 1) and <= Σ 1/2 1/2 - M2); see e.g., Ripley (1996 Section 2.1). We note that Rr is the smallest mis-classification error achievable by any classifier and.
8

---

Published in Transactions on Machine Learning Research (annuary/2000))
furthermore, is monotone decreasing as ||<|| increases. If || - 0 then RS - min 11 and Tp is no better than assigning every data point to the most prevalent class, while if ||<|| + co then Rr + 0 and Yr achieves perfect accuracy. The cases where Rp - 0 or Rs - min( 1 - are, theoretically, uninteresting and honce in this paper we only focus on the case where 0 < ||||| < co. We therefore make the following assumption.
Assumption 2. Let Se denote the set of indices i for which S * 0. Also let Co > 0, M > 0 and C₁ > 0 be constants not depending on such that 80 = |pg S M and
Ki 2 Co. JES
min
max (1x 1/2 pn |l ||-11/111 pall} S Cc
In addition, we assume that the number of spikes d' in the spiked covariance model (see Assumption ! fixed and does not grow with p or n. Remark 1. Assumption
is
g implies 0 Co S |ISI| S 2Cg- oc. We emphasize that sparsity is imposed on the whitened direction so that only few transformed features contribute to classification outcome; similar note that sparsity of the discriminanl direction == (1- (12) also implies that only a small subset of the raw covariates affects the response (as the classification boundary for a feature vector x is given by 1(x 3 > a} for some a E R). However, as noted by ZZhu & Bradic (2018) and Hall et al. (2014), spansity of 3 is incompatible with many real data application such as genome-wiee gene expression profiling where all genes are believed to play a role in discase markers, or analysis of micro-array data to identify leukemia or colon/pre pprostato cancer. In contrast sparsity of the whitened direction Ç = ( - 12) still allows for 3 = --1/2 be non-sparse, theoretical results and proofs, we assume that the sparsity level so is fixed and independent of the sample size n, though it may be arbitrary. The case where so grows with n can be handled mutatis mutandis. However, it requires more careful control of technical arguments and is therefore left to the interested reader.
assumptions can be found in Silin & Fan (022). To further explain this condition,
thereby circumventing the issue described above. Lastly, throughout
Our theoretical results are large-sample results in which the sample sizes nn and n2 for the training data increase the dimension p increases, and thus our next assumption specifies the asymplotic relationships between these quantities and the eigenvalues of >. Assumplion 3. Let be fixed and suppose that.
O(1), Inp o(n).
n2
Furthermore, the spiked eigenvalues , da of E satisfy
As = O(p),
for all & E [d].
If E satisfies Assumption 3 then tr(E) Op)) and r(E)- tr()/|| 0()). Recall that r(E) is the effective rank of Σ. Assumption
3 also implies that the spiked eigenvalues of Σ are unbounded as p increases
and this assumption distinguishes our theoretical results from existing results in the literature wherein it is
generally assumed that the eigenvalues of Σ are bounded; see Bickel & Levina (2004); Fan & Fan (2008) Shao
let al. 2011); b;; Cai & Zhang (209)) for a few examples of results under the bounded eigenvalues assumption. While the bounded eigenvalues assumption is prevalent, can also be problematic for high-dimensional data as it ignores the strong signals present many real data applications; see Fan et al. (2013) and Wang & Fan (017) for further discussions of this issue. Indeed, standard heuristic for PCA to keep the d largest principal components that explains 90% or 95% of the variability in the data, and hence if p is large and Σ has bounded eigenvalues then one has to choose d = (p) which is inconsistent with the use of PCA as a dimension reduction procedure.
If Σ has a spiked eigenvalue structure as in Assumption then W is given by
W= 1/2 4AA -II) o - uu)
(3:2)
and a natural estimate for W is
W = ú(A 022I 1/20 + 10--II--n -ûūT)
(3.3)
9

---

Published in Transactions on Machine Learning Research (January /2026
where Å is the diagonal matrix containing the d largest eigenvalues of the pooled covariance matrix, ú is the x d orthonormal matrix whose columns are the corresponding eigenvectirs, and 6² is as defined in Eq. (2.8). We now make the following assumption on U.
Assumption (Bounded Coherence). There is a constant Cu 2 independent of n and p such that
Cuvd
12
vP
Remark 2. The bounded coherence assumption appears frequently in statistical inference for matrix-valued data. More specifically, as every column 2/ has lz norm equal to 1, the rows U have, on average, & norm of order vd/p. Assumption 4 then guarantees that the maximum l norm of the rows U is also of order Vd/p. For more discussions about the bounded coherence assumption in the context of matrix completion, covariance matrix astimation, and random matrix theory, see Candès & Recht (2009): Fan et al. (2018): Rudelson & Vershynin (2016); Bloemendal et al. (2014) among others. Finally we note that 27 salisfies the bounded coherence assumption with constant Cu then
Cydrd p
Ea
+ +0²,
for all i E [p],
which logether with Assumption 3 implies max Ea O(1). In summary, Assumption allows for the spiked eigenvalues Ai.... Ad of the covariance matrix E to grow linearly with p while also guaranteeing that the entries of E remains bounded, i.e., each variable in Xis has finite variance.
We next state a result on the estimation accuracy of Ú. This result is a slight extension an earlier result in Cape et al. (2019). More specifically, Cape el al. (2019) assume E[X] = 0 and hence the sample covariance matrix simply " X X. In this paper we used the pooled sampled covariance matrix which requires first centering the feature vectors by the sample mcans of each class.
Theorem 1. Let X be a n p matric where the rows XI, are i.i.d samples from xNN(p1, E) + (1- ni (H2,, E) and E satisfies Assumption 1. Assumption & and Assumption 4 Let ů be the matriz of eigenvectors corresponding to the d largest eigenvalues of the pooled sample covariance matrix E. Then there crists d x d orthogonal matria Eu and constant C > O such that with probability at least 1 - O
d3 lp np
UEul22oo Sc
(3.4)
Note that, for simplicity, we assume in Theorem as well in the subsequent part of this paper that d, the number of spiked eigenvalues, is known. If is unknown then it can be consistently estimated using the ratio between consecutive eigenvalues, similar the procedures in Ahn & Horenstein (2013). More specifically, from Assumption 3 and Lemma we have with high probablitty that X = (p)) for S and Å O for 2 d+11 1. Here d are the eigenvalues of ε. There thus exists a significant. gap between d = 0(1) and Aa/Sa- = 2(p), and hence, letting be the smallest index k for which 2e/Ae41 2(1p)) we have d d= d asymptotically almost surely.
We now analyze the classification accuracy of lda o pea. Recall that the main idca behind Ida o pea is that we first construct an estimate Č for the whitened direction Ç and then project our whitened data onto the d largest coordinates, in magnitude, of Ć (see Algorithm [). By Eq. (3.1), the Bayes error rate Ry is a monotone decreasing function of ||||| and hence, to achieve Re is only necessary to recover the indices in S {i: & + 0}. In summary the error rate for lda pea converges to Rp provided we can (1) bound the error - { and (2) show that thresholding Č perfectly recovers Sç.
Theorem 2. Under Assumption 4 there exists a constant C > 0 such that with probability at least - O(ppp))
Inp n
1 -lloc S sC
(3.5)
10

---

Published in Transactions on Machine Learning Research (annuary/2200)
Assumption 2 implics &a * 0 if and only if S > Co for some constant Co O. By Theorem 2 if G 0 then Č = O( (n=11nn) with high probability. Let be hard thresholding of ñ, namely
- tn),
jEp
(3.6)
where in = (n-1 In p) for some 0 < a < 1.2 Given Ç, define the active set S = {j {3 03. The mis- classification rate for lda o pca. conditional on the training data {XI1,...,, Xim, } and {X11,.... Xim is then
Ridaopea = ((Iabel(Z) + T Idaopea (Z) {X Xims h. (X21, Xema
(3.7)
The following result shows that S recovers S exactly and Ida pea is asymptotically Bayes-optimal.
Theorem 3. Suppose that Z NTNNp(11, Σ) + (1 (122 E) where Tı (0, 1). Suppose Assumption Uhrough Assumption 4 are satisfied. We then have
P($ + Sc)= Opp 3)
(3.8)
Furthermore, suppose that lnp onn). We then have
Radaopen
Rpl 0
(3.9)
almost surely as n, ->
Our theoretical framework and results can be extended to elliptical distributions with minimal modifications; see the discussion in Section Key remarks our main results are presented below.
Remark 3. Theorem and Theorem 3 appear, at first blush, quite similar to the results in Mai et al. 2012 for
lassoed LDA using the discriminant direction 3 E- (11-1))))) -12).. However the underlying assumptions behind these results are very different. In particular lassoed LDA relies on irrepresentable condition analogous those in Zhao & Yu (2006) to achieve consistent variable selection, and thus the discriminative variables cannot be highly correlated with the remaining (irrelevant) variables. In contrast Theorem 3 do not assume an irrepresentable condition and furthermore Assumption !| only assumne boundedness of the non-spiked eigenvalues but allows for diverging spiked eigenvalues.
Remark 4. We now compare our results and settings with those in Sifaou et al. 2020a b. Firstly, Sifaou et al. (020aa) assumed: (1) proportional growth rate p/n c for some finite constant O, (2) distinct spiked eigenvalues Å > A2 > ... > Xa and lastly (3) ||ー|| = ((1). In contrast, our work is derived under a more challenging set-up wherein (1) the dimension and the number of training samples n satisfy n 'Ip + 0 (Assumption 3 (2) some or even all of the cigenvalues \ 2..
... ........ da > can be equal (Assumption I and (3) the leading eigenvalues { can diverge with n (Assumption allows for p= = n² for any >> 1, so that p/n -> 00.
1d
We emphasize that n -1 In +
The assumptions used in Sifaou et al. (020aa) because their thcoretical analysis rely heavily existing results in random matrix theory presented in Donoho & Ghorbani (2018); Donoho et al. 2018 which generally require p/n < 0o as well as bounded eigenvalues. In contrast, if p/n -> 00, then Donoho & Ghorbani (2018) shows that the best estimate of the true covariance matrix is typically a diagonal matrix (see also Bickel & Levina (2004)), which might make the results in Sifaou et al. (0000a) sub-optimal for the high-dimensional regime p/ n >> oo considered our paper.
Finally, Sifaou al. (2020a) did not guarantee that their classification rule is Bayes-optimal. Rather, they expression given in Theorem 3 of [Sifaou|
only show that the its missclaasiification rate converge to et al. (2020a), and is strictly larger than Bayes error when p/n increases. Siiaou et al. (2020b) subsequently extended the results in Sifaou et al. (2020a) to the case of heterogencous covariance matrices but the resulting classifier is once again not Bayes-oplimal. In contrast, Theorem 3 and Theorem 0 showed that the error rate for Ida pca and qda o pea (Quadratic Discriminant Analysis, QDA) converge to the Bayes error rate in the case of equal and unequal (class-conditional) covariance matrices, respectively. These results are, to the best of our knowledge, are among the first show Bayes consistency when combining PCA with LDA/QTI /DDA for high-dinensional classification under divergent spikes. A key reason for this improvement
11

---

Published in Transactions on Machine Learning Research (annuary/2000)
lies in two elements: (1) the assumption of sparsity in the signal vector Ç. which mitigates estimation error high-dimensional means and covariances (motivated by (Tibshirani et al. [2002]; Fan & Fan (2008)) and (2) the use of the entrywise matrix perturbation bounds, providing precise control over the estimated eigenspace ú and the resulting whitened direction Ç.
4 Numerical results
We now present simulation results and real data analysis for Ida o pea. Comparisons will be made against the nearest shrunken centroids method (NSC) of Tibshiraii al. 22002). sparse lincar discriminant analysis (SLDA) of (Shao et. al. (2011), direct sparse discriminant analysis (DSDA) of Mai et al. (2012), adaptive linear discriminant analysis (AdaLDA) Cai & Zhang (2019), and LDA rule with CAT scores (CAT-LDA)
of Zuber & Strimmer (009);; note that CAT-LDA uses different sphering transformation compared to Ida o pca. There are three other commonly used classifiers that are not included in our comparisons, namely the naive Bayes rule, the linear programming discriminant. (I.DD) (Cai & Liu 2011 and the regularized optimal affine discriminant (ROAD) (Fan et al. 2012). We omit these classifiers because (1) the naive Bayes rule special case of the NSC rule without soft thresholding of the mean vectors (2) Mai & Zou (2013)
showed that the ROAD and DSDA classifiers are equivalent and (3) Cai & Zhang (2019) showed that the
AdaLDA rule is a refinement of the LPD rule, i.e. compared to the LPD rule, the AdaLDA rule allows for "heteroscedastic constraints" and requires no tuning parameters.
Implementations of the DSDA and CAT-LDA rules are based on the TULIP and sda library in while the implementation of the AdaLDA rule is based on Matlab codes provided in ADAM github repository. We also consider the rotation pre-processing step of Hao et al. (2015). which yields the transformed data X11 N.....,4,,,, _ XIn, and X21, .aii....,,,,t X2n based on Matlab codes from HDRotation wheattttoonwwweeeeeeeeeee where
we set y = 0.25 in Eq. (2.10). These transformed data are then used as input to the DSDA and AdaLDA rules; we denote the resulting classifiers as DSDA(rot) and AdaLDA(rot).
4.1 Simulated examples
We consider three different simulation settings. For each setting, the number of features is set to p = 800 and we generate n2- = 100 data points from each class for the training data and also generate ni = 72 = 100 data points from each class for the testing data. The classification accuracy of the classifiers in each simulation setting are computed based on 200 Monte Carlo replications. The mean vectors of the two classes are
141 0s0o and μ22 (110: 0790); i.e., the vector μ2 contains 10 non-zero entries with values all equal to 1. We consider the following models for the covariance matrix Σ. These models were considered previously in Fan et al. 20122) Mai et al. 2011 Cai & Zhang 2019 Cai & Liu 2011), among others. Throughout this sub-section, the number of selected features refers to the dimensionality of the representation on which cach method operates. In particular, Ida o pca performs feature selection in the whitened feature space, whereas other methods select features directly from the raw feature space. See Section 7.2 for a detailed study in which the sparsity levels of the mean difference μ22 - H1, the discriminant direction 3, and the whitened direction Ç are matched.
Model 1 (equal correlation): Here Σ= = (ax))xx = p11 + (1- P) i.e., Oij = for i * j and Oij = 1 for i i= j. With this covariance structure Σ, the discriminant direction / E E 14) and whitened direction (= E-=- - no) are non-sparse (all entries of and Ç are non-zero).
2. Model 2 (block diagonal with equal corvelation) Here Σ = (ar)) is assumed to be a block diagonal matrix with two blocks of size 20 x 20 and (p 20) x (p 20). Both diagonal blocks are also of the form p111 + (1-p)I where the correlation p the same for both blocks. The discriminant direction and whitened direction Ç are sparse in this model, with B3 and ( both having 20 non-zero entries.
Model 3 (random correlation) Here (O 1) and ez = minteji [CC Note that we generate a new L for every Monte Carlo replicate. For further comparison, we also consider Lij generated from the uniform distribution on [-1,1]
LLT + coTP where E Rpx10 with Luj generated from
12

---

Published in Transactions on Machine Learning Research (January /2026
Table 1: Mis-classification rate (%) with standard deviations (%) in parentheses for equal correlation setting
(model 1), based on 200 independent Monte Carlo replicates.
0.50
0.60
0.70
0.80
0.90
Oracle
1.37 (0.87)
0.61 (0.53)
0.19 (0.30)
0.19 (0.30)
0.00 (0.00)
Ida pca.
1.74 (1.00)
1.00 (0.82)
0.55 (0.67)
0.55 (0.67)
0.22 (0.39)
CAT-LDA
7.18 (2.22)
4.70 (1.77)
2.60 (1.34)
1.10 (0.88)
0.38 (0.56)
DSDA
3.27 (1.47)
1.96 (1.08)
0.74 (0.62)
0.76 (0.65)
0.00 (0.04)
DSDA
6.31 (1.83)
3.26 (1.27)
1.16 (0.83)
0.30 (0.41)
0.01 (0.06)
AdaLDA
4.15 (1.61)
2.62 (1.24)
1.22 (0.91)
0.31 (0.42)
0.00 (0.04)
(rat)
7.04 (2.07)
3.70 (1.26)
1.19 (0.79)
0.13 (0.26)
0.00 (0.00)
SLDA
17.97 (3.53)
14.46 (2.99)
11.14 (2.43)
11.14 (2.43)
1.64 (0.92)
NSC
20.38 (8.53)
22.60 (8.38)
24.64 (8.36)
24.64 (8.36)
29.26 (8.04)
and the Student's t-distribution with 5 degrees of freedom. The discriminant direction 8 and Ç are generally non-sparse in this model, i.e. all of their entries are non-zero.
The tuning parameters for each classifier are chosen using five-fold cross validation (CV). In particular, SLDA (Shao et al. 2011) requires two tuning parameters, one being the number of non-zero entries in fa - /12 and the other being the number of non-zero entries in E. For simplicity we shall assume that the sparsity of AI 12 is known when implementing SLDA and thus the only tuning parameter required is the number of non-zero entries in E; following Cai & Zhang (2019), this tuning parameter is selected from the set of values {vn--lu 1.5nn Inp...., 5 n- 1 In p} using five-fold CV. The Ida pca classifier in Algorithm [ also requires two tuning parameters, namely (1) the number of spikes d in the estimation of Σ and (2) the sparsity level 8 in Ç. We chose d to account for at least 90% of the total variability in the data, i.e., is the smallest integer of satisfying (Etz {))tr(E) 2 0.9; hore a 2 du 2 ... are the cigenvalues of the pooled sample covariance matrix E. We acknowledge that the choice of d can significantly impact performance; therefore, a sensitivity analysis is provided in Section 7.2 For each Monte Carlo replication, the training data are randomly split into five folds for CV, with new splits generated each time. The sparsity level s is chosen from s E f1.2...., 30) to minimize the average CV misclassification error. In the case where multiple values of s yield the same minimum error, we choose the smallest such 8 to promote a more parsimonious model. Similarly, the number of features in CAT-LDA selected among the top {1, ,, features ranked by CAT scores using five-fold CV. The choice of the upper limit 30 is motivated by the common sparsity condition 80 ((nn/ In p)",, for some 7 > O, for instance, T 1/2 in Cai & Liu (201T). We deliberately choose n/ In (= 30) as a relaxed upper bound to broaden the grid and ensure adequate coverage.
For data generated according to Model 1, Table [ and Table 2 show that Ida o pca achieves the highest accuracy while using only a small number of features compared to the other classifiers. Note that although CAT-LDA also apply a whitening transformation before performing LDA, the accuracy of CAT-LDA is much worse compared to that of lda pca. The DSDA and DSDA classifiers have slightly better accuracies compared to those for the AdaLDA and (roL classifiers; recall that the (SDA(rot)) and AdaLDA(rot)) classifiers first applied the rotation pre-processing step of Hao et al. (2015) before running DSDA and AdaLDA on the transformed data. The NSC classifier has the largest cissiiissiiinntion error; this is a consequence of the NSC rule ignoring the correlation structure in Σ. Finally, the oracle classifier correspond to the LDA rule where IT, Pl, μ1 and Σ are known, and thus its error rate is the Bayes error rate from Eq. (.11).
Remark 5. Rocall that for Model 1, all entries of the discriminant direction B and the whitened direction Ç are non-zero. These entries however can be classified into those representing strong signals vs weak signals based on their magnitudes; see Table the supplementary. The strong signals appear in the first 10 entries of B (similarly <) and the remaining p- - entries of B (similarly () correspond the weak signals. For example, if p p= 0.5 then the first 10 entries of 3 are all equal 1.98 and the remaining p p= 10 entries are all equal to -0.02. From Table 2 we see that. lda pea kept all of the coordinates corresponding to the strong
13

---

Published in Transactions on Machine Learning Research (January /2026
Table 2: Average number of nonzero coefficients with standard deviations in parentheses for equal correlation sctting (model 1), based 200 independent Monte Carlo replicates.
0.50
0.60
0.70
0.80
0.90
Ida pea
12.04 (4.53)
11.31 (4.10)
9.52 (4.01)
9.52 (4.01)
3.68 (0.98)
CAT-LDA
24.85 (3.74)
24.70 (3.58)
24.60 (3.97)
22.29 (3.93)
15.45 (2.39)
DSDA
96.20 (31.06)
106.42 (32.57)
117.15 (30.52)
117.05 (29.97)
96.48 (8.66)
DSDA(rot)
33.78 (32.71)
36.88 (38.13)
57.23 (61.02)
147.57 (66.97)
176.29 (5.08)
AdaLDA
46.54 (5.95)
45.84 (3.25)
46.72 (6.41)
47.86 (5.88)
48.65 (5.02)
(rot)
5.68 (1.82)
5.57 (1.84)
5.86 (1.82)
7.83 (2.54)
19.85 (9.96)
SLDA
728.84 (226.63)
788.16 (95.69)
799.96 (0.21)
799.96 (0.21)
799.96 (0.21)
signals and only added a few coordinates corresponding to the weak signals; this explains the small number of features used in Ida pea. In contrast, the DSDA and SLDA classifiers include a large number of (noisy) features with weak signals.
For Model 2, Table 3 shows that SLDA performs slightly better than Ida pea and DSDA. However, from Table we see that SLDA also selects much larger number of features compared to both Idaopca and DSDA. Recall that, for model 2, both 3 and ( contains exactly 20 non-zero entries and hence SLDA is selecting large number of extrancous, non-informative features; a similar, albeit much less severe, phenomenon is observed for DSDA. Table GaableBare and shows that Idaopca, CAT-LDA, and AdaLDA have comparable accuracy with a similar number of selected features. The NSC rule once again has the largest mis-classification error due to it ignoring the correlation structure in Σ. Finally we see that the rotation pre-processing step described in Hao et al. (2015) lead to substantial loss in accuracy for the DSDA and AdaLDA classifiers; to understand why this happens, we need to extend the theoretical analysis in Hao et al. 2015 (which assume that Σ and Mpot.p are known) to the setting where E and Urorpp have to be estimated.
We leave this investigation for
future work.
Finally for Model 3, Table 5 and Table G show that Ida pea has both the highest accuracy as well as the smallest number of selected features among all the considered classifiers; DSDA has a slightly worse accuracy and also selected much larger number of features, when compared to Ida pca. SLDA now has the worst accuracy and also selects almost all p 800 features, and this is a consequence of Σ being a dense matrix. The mussclassification rate for AdaLDA and AratDA(rot) are also quite large, and surmise that this is due to the numerical instability when solving the linear programming problem in AdaLDA. Indeed, the condition numbers for can be quite large; see Table 12 in the supplementary for summary statistics of these condition numbers using the same 200 Monte Carlo replicates as that for generating Table 5 Following the suggestion in Cai & Zhang (209)) we replace the sample covariance matrix E used in the optimization problem for AdaLDA and AdaLDA(root (see Eq. (8) and Eq. (9) in Cai & Zhang (2019)) with E = E + vn- Inpppnn Ip: the resulting classifiers are denoted as AdaLDA(reg) and AdaLDA(rot + reg), respectively. Table 5 however shows that using Σ only leads to minimal increase in accuracy. Finally we observe that the rotation pre-processing step once again leads a substantial loss in accuracy for the DSDA and AdaLDA classifiers.
Remark 6. While our ical results primarily focus on fixed sparsity (see Remark [ for more growing sparsity), we deliberately include dense examples with so p, such as Model 1 and Model 3, to demonstrate robustness beyond the sparse regime. In the dense models, full recovery of Ç or 3 is not anticipated but Ida o pea still delivers strong classification performance due to ellective decorrelation. In particular, for Model 1, our method's performance is very close to the oracle across varying p.
4.2 Real data analysis
We now assoss the performance of lda o pea on two gene expression data sets for leukaemia and lung cancer. The leukaemia dataset (Golub et al. 1999) includes p = 7128 gene measurements on 72 patients with either acute lymphoblastic leukemia (ALL)) acute myeloid leukemia (AMT.) and wish to classify patients into
14

---

Published in Transactions on Machine Learning Research (January /2026
Table 3: Mis-classification rate (%) with standard deviations (%) parentheses for block diagonal setting
(model 2), based on 200 independent Monte Carlo replicates.
0.50
0.60
0.70
0.80
0.90
Oracle
5.30 (1.71)
5.21 (1.52)
5.34 (1.64)
5.48 (1.62)
5.40 (1.74)
Ida pca.
9.08 (2.71)
8.61 (2.59)
8.30 (2.55)
8.26 (2.45)
8.12 (2.58)
CAT-LDA
10.08 (3.00)
9.36 (2.57)
9.31 (2.50)
9.25 (2.59)
9.21 (2.89)
DSDA
9.63 (2.61)
9.62 (2.64)
9.50 (2.66)
9.58 (2.67)
9.47 (2.42)
DSDA
20.64 (4.16)
20.52 (3.99)
20.75 (4.18)
21.02 (4.44)
20.90 (4.40)
AdaLDA
12.70 (2.99)
13.06 (3.39)
13.88 (3.45)
13.75 (3.48)
14.68 (3.86)
AdaLDA(rot)
23.91 (3.79)
23.81 (3.72)
24.33 (4.10)
23.89 (3.92)
24.12 (3.94)
SLDA
6.99 (2.89)
6.66 (1.92)
6.71 (2.17)
6.62 (1.97)
6.56 (1.90)
NSC
25.14 (3.11)
26.47 (3.31)
28.61 (3.35)
29.59 (3.24)
30.52 (3.82)
Table 4: Average number of nonzero coefficients with standard deviations parentheses for block diagonal setting (model 2), based on 200 independent Monte Carlo replicates.
0.50
0.60
0.70
0.80
0.90
Ida pea
20.19 (4.87)
20.48 (5.06)
20.48 (4.88)
20.96 (4.96)
20.58 (4.44)
CAT-LDA
20.02 (4.95)
19.86 (4.31)
20.05 (4.62)
20.32 (4.21)
20.33 (4.54)
DSDA
50.28 (26.61)
50.74 (28.92)
47.48 (25.77)
43.66 (22.42)
45.15 (23.78)
(rDA
51.02 (39.23)
46.13 (29.09)
51.93 (36.08)
45.56 (31.15)
49.62 (33.50)
AdaLDA
19.45 (2.37)
18.53 (2.23)
17.93 (2.27)
17.40 (2.20)
16.85 (2.07)
(rot
14.40 (3.14)
14.81 (3.53)
15.16 (3.40)
15.69 (3.99)
14.97 (3.84)
SLDA
405.98 (390.84)
320.17 (380.33)
210.99 (336.12)
242.19 (352.81)
172.00 (309.62)
either ALL or AML based their gene expressions. The training set consists of the gene expressions for Ri 27 patients with ALL and n2 = 11 patients with AML while the test set contains the gene expression data for 20 patients with ALL and 14 patients with AML. The lung cancer data was originally analyzed in Gordon et al. (2002) and we use version of the data wherein the predictor variables with low variances are removed; see Pun & Hadimaja (2021). The resulting data set contains tumor tissues with = 1577 features collected from patients with adenocarcinoma (AD) or malignant pleural mesothelioma (MPM). According Gordon al. (2002), MPM is highly lethal but distinguishing between MPM and AD quite challenging in both clinical and pathological settings. The training set consists of nı 120 gene expressions for patients with AD and n2 25 gene expressions for patients with MPM, and the test set consists of gene expressions for 30 patients with AD and patients with MPM. Table and Table X presents the classification accuracy and number of selected features for various classifiers when applied to the leukemia and lung cancer data, respectively. The hyperparameters for these classifiers are selected using leave-one-out CV. Table 7 and g indicate that the performance Ida o pea competitive with existing state-of-the-art classifiers while also operating on substantially lower siverddmensiooall representation.
5
Extensions of Ida o pca
5.1 Multi-class classification
Suppose given training data from K > 3 classes where the feature vectors for each class are iid samples from p-dimensional multivariate normal distribution, i.e.,
XuI...,, Xin, No(ui 2), i [K].
15

---

Published in Transactions on Machine Learning Research (January /2026
Table 5: Mis-classification rate (%) with standard deviations (%) in parentheses for random correlation setting (model 3), based on 200 independent Monte Carlo replicates.
Method
u(-1,1)
N(0,
ls
Ida pca
5.07 (2.40)
12.39 (4.17)
13.72 (5.00)
CAT-LDA
11.93 (4.04)
23.67 (6.13)
25.68 (6.82)
DSDA
7.32 (3.05)
16.92 (5.21)
17.90 (5.49)
DSDA(rot)
18.51 (6.37)
31.81 (7.11)
33.65 (7.28)
AdaLDA
17.94 (6.03)
44.55 (7.12)
47.89 (5.10)
(rot.
23.00 (7.49)
39.29 (6.11)
42.52 (6.28)
AdaLDA(reg)
29.48 (7.58)
47.49 (4.52)
48.96 (4.04)
AdaLDA(reg
21.23 (7.43)
38.93 (6.60)
42.52 (6.28)
SLDA
50.10 (4.40)
49.98 (4.25)
50.09 (3.70)
NSC
30.30 (6.28)
44.55 (5.16)
46.70 (4.49)
Table 6: Model size with standard deviations in parentheses for random correlation setting (model 3), based 200 independent Monte Carlo replicates.
Method
u(-1,1)
N(O, 1)
t.5
Ida o pca CAT-LDA
11.93 (4.10)
11.48 (3.88)
11.37 (3.92)
21.19 (6.51)
20.64 (6.88)
20.75 (6.90)
DSDA
66.00 (31.92)
65.63 (34.45)
64.91 (30.51)
DSDA(rot)
51.08 (34.82)
71.94 (41.64)
77.56 (49.29)
AdaLDA
15.99 (2.32)
10.46 (1.97)
8.51 (2.22)
AdaLDA(rot
14.62 (4.61)
19.08 (3.93)
16.26 (3.43)
AdaLDA (reg)
11.38 (1.77)
8.50 (2.00)
7.12 (1.96)
(regDA(reg-rrof)
9.86 (4.71)
15.99 (3.63)
16.26 (3.43)
SLDA
799.51 (1.31)
800.00 (0.07)
800.00 (0.00)
NSC
12.21 (8.47)
55.38 (131.65)
102.89 (192.83)
Here n, denote the number of training data points from class E [K]. The testing data point Z is drawn from a mixture of K multivariate normal distributions namely, Z Σ m,Np(u,; Σ) with क, 2 and क. Σm 1. Note that, for case of exposition, the (Xij} are assumed to be multivariate normals but the subsequent results also hold when } are elliptically distributed as in Section 7.3
Now define, for 2 < i S K, the whitened direction ((i) and the whitened indices S, via
( = -W(ss m1),
S; (j 5500 #0}.
Define the global whitened set as S. S₂ U. .... U SK. The indices in Se are the important variables for feature selection. The extension of Ida pea K classes is then given in Algorithm below, and theoretical results presented in the supplementary.
16

---

Published in Transactions on Machine Learning Research (January /202)
Table Misccaassiiicttion rate and model size of various methods for the leukaemia data
Method
Training Error
Test Error
Model Size
Ida o pca
0/38
1/34
12
DSDA
0/38
1/34
36
AdaLDA (reg)
0/38
1/34
18
NSC
1/38
3/34
24
Table 8: Mis-classification rate and model size of various methods for the lung cancer data
Method
Training Error
Test Error
Model Size
Ida o pca.
0/145
0/36
20
DSDA
1/145
0/36
25
AdaLDA (reg)
0/145
0/36
388
NSC
1/145
0/36
1206
Algorithm 3: K-classes Ida pca decision rule
Input: X,, i [K], E and the test sample Z Output: Idaopeca
Algorithm
// Step 1: Perform PCA on the feature vectors for the training data (standard PCA approach).
Extract the d largest principal components obtain W.
// Step 2: For 2 i < K, let X; WX, and ((0) = X; - X,. Also let S be the set of indices corresponding to the 8; largest elements of g(0) in absolute values. The value of 8; is, general, a user-specified or tuning parameter. Nevertheless, under certain conditions, we can also estimate sç see Eq. (855] in the supplementary.
201) the vector obtained from // Step 3: Given the test data point Z, let Ž WZ and denote by Ss
(00 by keeping only those coordinates belonging to S: // Step 4: Set D = and calculate, for 2 < S K, the discriminant score for class relative to class 1 as
D. = 2 -
x
S, Sg.
201 In
+
74
(5.1)
n2
// Step 5: Assign the label of to the class that maximizes the discriminant score, i..,,
T Idaopca (Z) arg max D..
(5.2)
.E[K]
5.2 Heterogeneous covariance matrices
Despite the simplicity and popularity of regularizod or sparse LDA for high-dimensional data, the assumption of equal covariances not always tenable in practice. More specifically, suppose we are given a p-variate random vector Z drawn from mixture TNN Σn) + W E2) with En possibly distinct from
17

---

Published in Transactions on Machine Learning Research (January /2026)
Σ.. The Bayes classifier is then the QDA rule given by
[1 if (Z) (z p1) -(z p2) 'I2' (Z p2) S (5.3)
Tz
2 otherwise
where x = 21n T1 In ΣЛ = Σ2, then Eq. 153 reduces to Eq- 1.1 We now discuss how the I-
rosults in Section can be extended to quadratic discriminant analysis (QDA) with PCA. Firstly assume that En and Σ② both have spiked covariance structures specified below.
Assumplion 5. Let uil,...; Wid,, for = 1, 2 be orthonormal vectors in IR" and assume that the covariance matrix Σ, for the p-variate normal distributions Ny Σx) is the form
Σ;
Au uuu i +o =UAU + +o 1=
(5.4)
Here A; diag(Au: d; d: diagonal matrix and U4 p x d; matrix with orthonormal columns. We assume implicitly that Xa 2 ... 2 Aid 0, > 0 and d, < p, = 1,2.
Note that a recent line of research on QDA for high-dimensional classification is based on the assumption
that Σn - -E- sparse; see e.g., Li & Shao (2015); Jiang et al. (2018); Cai & Zhang (021)). In contrast, Assumption do not enforce any sparsity assumption and also allows for ΣS and Σg to have different spiked cigenvalues and cigenvectors; the latter is generalization of the common principal components (CPC)
assumption in Zhu 20060; Pepler et. al. (2017): Flury 1988) where the leading principal components the
same for both 11 and 2g.
Under Assumption the whitening transformation W: = 1/2 for i= 1,2 is of the form in Eq. 2.7 and thus a suitable estimate for W; is given by Eq. (3.3) Let C = W,y for i= 1,2, and denote the whitened index sets by A4 = [j Sij * 0} for i = 1,2. Let Ag = Ai U A2 and note that the elements in Ao are the signal coordinates for the QDA rule (after the PCA step). Eq. (5.3) can be written as
T(Z) = 1
if [[M (Z- -11))aa1011 1l|M2 #2)1 12 12 2 otherwise
(5.5)
A plugin decision rule is then obtained by replacing μi, Eg and Ao with their estimates X,, Sr and Ao Note that the intercept is non-trivial to estimate in the high-dimensional setting as it involves the log- determinant In E. see Cai et al. (2015) for further details. For our paper employ the data-driven approach of Jiang et al. (2018) which circumvents the need to estimate the determinants of Σ, and Σg. The full details of qda o pea are specified in Algorithm below, and theoretical results for qda o pea are provided in the supplementary.
18

---

Published in Transactions on Machine Learning Research (January /2026
Algorithm 4: qda pea decision rule
Input: X1, X2, Σi, Eg and test data point Z
Output: f Tddacpoa
Algorithm
// Step 1: Perform PCA the feature vectors for the training data. (standard PCA approach). Extract, for i = 1,2, the dj largest principal components UG and compute Wi as in Eq. (3.3).
// Step 2: For = 2, set & = W;X; and form the indices set A by selecting the st largest (in magnitude) coordinates of G,
// Step Le As = -r U Åg and define for any æ E R2
Q(x X2,
X2, W,, W2 Aod) | x1) 112 x2)
// Step Finind & to minimize the empirical (-1 loss of the decision rule induced by Q, i.c.,
1
arg min WER NE
((T(X)) X2, M. We, Mhhresh)) ++)+ i)
(5.6)
where, for any Vhhresh E ER, we define
[ if Q(1y) (X1, otherwise
S Thhresh
M(X,, (81, X2, Wi, W2, Ao, (hhresh))
// Step 5: Given a test data point Z, return the decision rule
1 if Q(Z) W,, W2 Aot) S R otherwise
qdaopea (Z)
(5.7)
As indicated by Jiang et al. (0188). Thhresh is selected by minimizing the in-sample misclassification error based on X1, X2, Σ1, and Σ2. Under Assumption 5 the grid search can be narrowed to a neighborhood around log |i21
This is justified by the improved stability of the log-determinant ratio in the spiked setting, where bulk eigenvalue of is consistently estimated via the pooled sample covariance and the leading spikes remain well-separated (from the bulks).
log A1
7iz
6 Discussion
In this paper we addressed the classification problem for high-dimensional data by analyzing the prototypical Ida pca classifier that first transforms the feature vectors using a whitening transformation, then performs feature selection on the whitened data, and finally applies LDA in the dimensionally reduced space. We show that, under spiked eigenvalue structure for Σ, the mis- dus-classification error rate for Ida pea is asymptotically Bayes optimal whenever + co and n Inp + O. While the Bayes consistency of lda pca is similar to that of classifiers based on estimating the discriminant direction 3, the underlying assumptions and motivations for our results are substantially different. Indeed, the focus on PCA and the whitening matrix leads to the natural assumption that Σ has spiked/diverginn /iverging cigenvalues while earlier results that focus on estimation of 3 had generally assumed that E is sparse or that the eigenvalues of Σ are bounded. Numerical experiments indicate that Idaopea is competitive with existing slate-of-the-art classifiers while operating substantially
19

---

Published in Transactions on Machine Learning Research (January 2202)
lower-dimensionll representation. This behavior persists even when the underlying sparsity levels in the raw and whitened feature spaces are matched; see more details in Section 7.2
We now mention two interesting issues for future research. The first is to extend the theoretical results in this paper for combining LDA with PCA to other, possibly non-linear, dimension reduction techniques such as (classical) multiddimensional scaling, kernel PCA, and Laplacian eigenmaps (Belkin & Niyogi 2003), followed by learning classifier in the dimensionally reduced space, The second issue concerns the spiked covariance structure in Assumption [ In particular, while Assumption [1] is widely used in the literature, see c.g., HHo the non-spiked cigenvalues and eigenvectors might be somewhat restrictive. We can consider relaxing Assumption by assumiin that Σ arises from an approximate factor model (Fan et al. 2013) or that Σ can be decomposed into a low-rank plus sparse matrix structure (Agarwal et al.) 2011). We surmise that, due to the focus on the whitening matrix Σ-1/2, theoretical analysis of Ida pea under these more general covariance structure will also leads interesting technical developments; c.g., while perturbation results for Đ-1 E-- given E E are well-studied, much less is known about perturbation bounds for S-1/2 E-1/2 given Σ - E.
et al. 2015;; Cai et al. (2013); Johnstone (2001): Birnbaum al. (2012). its assumption
References
Alekh Agarwal, Sahand Negahban, and Martin Wainwright. Noisy matrix decomposition via convex re- laxation: Optimal rates in high dimensions. Annals of Statisties, 40:1171-1197, 2011. doi: 10.1214/ 12-AOS1000.
S. C. Ahn and A. R. Horenstein. Eigenvalue ratio lest for the number of factors. Econometrica, 81:1203 1227, 2013.
Lucia Alessi, Matteo Barigozzi, and Marco Capasso. Improved penalization for determining the number of factors in approximate factor models. Statistics Probability Letters, 80(23-24):1806 1813, 2010.
Theodore Wilbur Anderson. An Introduction to Multivariate Statistical Analysis, Third Edition. New Jersey Intersecnt,,,,,ecee, 2003. ISBN 9780471360919.
Jushan Bai and Serena Ng. Determining the number of factors in approximate factor models. Econometricn, 70(1):191-22) 2002.
P. N. Belhumeur, J. P. Hespanha, and D. J. Kriegman. Eigenfaces vs fisherfaces: recognition using class specifie linear projection. IEEE Tronsactions on Pattern Analysis and Machine Intelligence, 19:711 720, 1997.
M. Belkin and P. Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation, 15:1373-133 1396, 2003.
Quentin Berthet and Philippe Rigollet. Optimal detection of sparse principal components high dimension. Annals of Statistics, 41:1780 1815, 2012. doi: 11221413-AAS111
Peter Bickel and Elizaveta Levina. Some theory for Fisher's linear discriminant function, 'naive Bayes', and some alternatives when there are many more variables than observations. Bernoulli, 10:989-1010, 2004.
Peter Bickel and Elizaveta Levina. Covariance regularization by thresholding. Annals of Statistics, 36: 2577-2604, 02 2009. doi: 400555000-000000
AOS600.
Xin Bing and Marten Wegkamp. Optimal discriminant. analysis in high-dimensional latent factor models. Annals of Statistics, 51(3):1232 1257, 2023.
Xin Bing, Bingging Li, and Marten Wegkamp. Linear discriminant regularized regression. arXio preprint arXiv:2402.14.00 2024.
Aharon Birnbaum, lain Johnstone, Boaz Nadler, and Debashis Paul. Minimax bounds for sparse PCA with noisy high-dimensional data. Annals of Statisties, 41:1055-1088 1084, 2012. doi: 10.1214/12-A0S1014
20

---

Published in Transactions on Machine Learning Research (January /2026
A. Bloemendal, L. Erdös, A. Knowles, H.-T. Yau, and J. Yin. Isotropic local laws for sample covariance and generalized wigner matrices. Electronic Journal of Probabilily, 19:1-53, 2014.
T. Cai and Weidong Liu. A direct estimation approach to sparse linear discriminant analysis. Journal of the American Statistical Association, 106:1566 -5777, 2011. doi: 20.1199/jaa.2211..mm...
T. Cai and Linjum Zhang. High-dimensional linear discriminant analysis: Optimality, adaptive algorithm, and missing data. Journal of the Royal Statistical Society, Series B., 81:675-705, 2019. doi: 10.1111/rssb. 12326.
T. Cai, Zongming Ma, and Yihong Wu. Optimal estimation and rank detection for sparse spiked covariance matrices. Probability Theory and Related Fields, (61:781-815) 2013. doi: 10.1007/00044000144400022..
T Tony Cai and Anru Zhang. Rate-optimal perturbation bounds for singular subspaces with applications high-dimensional statistics. Annals of Statistics, 46(1):60 89, 2018.
T Tony Cai and Linjum Zhang, A convex optimization approach to high-dimensional sparse quadratic discriminant analysis. Annals of Statistics, 49(3):1537-1568, 2021.
T Tony Cai, Tengyuan Liang, and Harrison Zhou. Law of log determinant of sample covariance matrix and optimal estimation differential entropy for high-dimensional gaussian distributions. Journal of Multivoriate Analysis, 137:161 172, 2015.
Emmanuel Candes and Terence Tao. The Dantzig selector: Statistical estimation when is much larger than n. Annals of Statislics, 35:2313 2351, 2005. doi: 10.1214/0090333600000000
Emmanuel Candès and Benjamin Recht. Exact matrix completion via convex optimization. Communications of the ACM, 9:717-772, 11 2009. doi: 10.1007//112008-009 904--5.
Joshua Cape, Minh Tang, and Carey Priebe. The two-to-infinity norm and singular subspace geometry with applications to high-dimensional statistics. Annals of Statistics, 47:2405- 24399, 2019. doi: 10.1214/ A0S17555
David Donoho and Behrooz Ghorbani. Optimal covariance estimation for condition number loss in the spiked model. arXiv preprint orXiv:1810.074003, 2018
David Donoho, Matan Gavish, and lain M Johnstone. Optimal shrinkage of eigenvalues in the spiked covariance model. Annals of statistics, 46(4):1742, 2018.
Junbo Duan, Charles Soussen, David Brie, Jérôme Idier, Mingxi Wan, and Yu-Ping Wang. Generalized lasso with under-determined regularization matrices. Signal processing, 127:239 246, 2016.
Y. C. Eldar and A. Oppenheim. MMSE whitening and subspace whitening. IEEE Transactions on Information Theory, 49(7):1846-185) 2003. doi: 10.1109/TII 201TT22003813500.
Jianqing Fan and Yingying Fan. High dimensional classification using features annealed independence rules. Annals of Statistics, 36:2605-2637 2637, 02 2008. doi: 10.1214/07-A0S504.
Jianqing Fan, Yang Feng, and Xin Tong. A road to classification in high dimensional space. Journal of the Royal Statistical Sociely, Series B, 74:745 771, 09 2012. doi: 10.1111/j.1467-990 2012.01029.x
Jianqing Fan, Yuan Liao, and Martina Mincheva. Large covariance estimation by thresholding principal orthogonal complements. Journal of the Royal Statistionl Sociely, Series B, 75:603 680, 2013. doi: 10.
1111/rssb.12116
Jianqing Fan, Weichen Wang, and Yiqiao Zhong, An Eoo eigenvector perturbation bound and its application to robust covariance estimation. Journal of Machine Learning Research, 18:1-42, 03 2018.
Jianqing Fan, Kaizheng Wang, Yiqiao Zhong, and Ziwei Zhu. Robust high dimensional factor models with applications to statistical machine learning. Statistical Science, 36(2):303, 2021.
21

---

Published in Transactions on Machine Learning Research (January /2026
Kai-Tai Fang and T.W. Anderson. Statistical inference in elliptically contoured and related distributions.
Journal of the Americon Statistical Associalion, 87:253-2555 254, 03 1992. doi: 10.2307/2290488
Xu Feng and Wenjian Yu. A last adaplive randomized pea algorithm. In IJCAI, pp. 3695 3704, 2023.
Bernhard Flury. Common Principal Components & Related Multivariate Models. John Wiley & Sous, Inc., 1988.
T.R. Golub et al. Molecular classification of cancer: class discovery and class prediction by gene monitoring, Science, 286:531-7, 11 1999.
Gavin Gordon et al. Translation of microarray data inlo clinically relevant cancer diagnostic tests using gene expression ratios in lung cancer and mesothelioma. Cancer Research, 62:4963-1 7, 10 2002.
Nathan Halko, Per-Gunnat Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. STAM review, 53(2):117-288, 2011.
Peter Hall, Jiashun Jin, and Hugh Miller. Feature selection when there are many influential features.
Bernoulli, 20(3):1647 1671, 2014. doi: 10.3150/13-BEJ5366
Marc Hallin and Roman Liška. Determining the number of factors in the general dynamic factor model. Journal of the Americon Statistical Association, 102(478):603 617, 2007.
Ning Hao, Bin Dong, and Jianqing Fan. Sparsifying the Fisher linear discriminant by rotation. Journal of the Royal Statistical Society, Series B, pp. 827 851, 2015. doi: 10.1111/rssb.120 12092.
R. A. Horn and C. R. Johnson. Topics in Matrix Analysis. Cambridge University Press, 1991.
Binyan Jiang, Xiangyu Wang, and Chenlei Leng. A direct approach for sparse quadratic discriminant analysis. Journal of Machine Learning Research, 19:31:1-31 31:37, 2018.
lain Johnstone. On the distribution of the largest eigenvalue in principal components analysis. Annals of Statistics, 29:295 327, 2001. doi: 10.1214/a0s/ aa////00//11000////
lan Jolliffe. Principal Component Analysis. Springer Verlag, New York, 2002.
lan T Jolliffe. A note on the use of principal components in regression. Journal of the Royal Statistical Society, Series C., (3):3000 303, 1982.
Thibaut Jombart, Sébastien Devillard, and François Balloux. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC genetics, 11(1):1-15, 2010.
Agnan Kessy, Alex Lewin, and Korbinian Strimmer. Optimal whitening and decorrelation. The American Statistician, 72:309-311 314, 2015. doi: 00.1000000031105.200000 2010:311135.2...11211 2016.1277159
Vladimir Koltchinskii and Karim Lounici. Concentration inequalities and moment bounds for sample covari- ance operators. Bernoulli, 23:110 133, 2017. doi: 1..1150/15-BEJ7300
Quefeng Li and Jun Shao. Sparse quadratic discriminant analysis for high dimensional data. Statistica Sinica, pp. 457 473, 2015.
Karim Lounici. High-dimensional covariance matrix estimation with missing observations. Bernoulli, 20: 1029-1058, 2014. doi: 1-3150/12-BBJ453
Qing Mai and Hui Zou. A note on the connection and equivalence of three sparse linear discriminant analysis methods. Technometrics, 55:243-2444 4246 2013. doi: 10.1080 10.1080/004017000 700.006...76600
Qing Mai, Hui Zou, and Ming Yuan. A direct approach to sparse discriminant analysis in ultra-high dimen- sions. Biometrika, 99:29 42, 2012. doi: 10.2307/41720671
Qing Mai, Yi Yang, and Hui Zou. Multiclass sparse discriminant analysis. Statistica Sinica, 29:97-111, 2015. doi: 10.5705/ss. 2020166.11.
22

---

Published in Transactions on Machine Learning Research (annary//020)
Yue Niu, Ning Hao, and Bin Dong. A new reduced-rank linear discriminant analysis method and its appli-
cations. Stalistion Sinica, 28:189-28 202,2002-22000 2015. doi: 10.37055ss.2222155.00 5.0387.
P Theo Pepler, DW Uys, and DG Nel. Discriminant analysis under the common principal components and Computation, 4666):4812-48277, 2017.
model. Communications in
M. M. Prasad, M. Sukumar, and A. G. Ramakrishnan. Orthogonal LDA in PCA transformed subspace.
In Twelfth International Conference on Frontiers in Handwriting Recognition, pp. 172 175, 11 2010. doi: 10.1 WW009//CFH.2000000
Chi Seng Pun and Matthew Hadimaja. A self-calibrated direct approach to precision matrix estimation and linear discriminant analysis in high dimensions. Computational Statistics and Data Analysis, 155:107105, 03 2021. doi: 10.1016/j.eee 2at...2220.1100000
B. D. Ripley. Pattern Recognition and Neural Networks. Cambridge University Press, 1996.
M. Rudelson and R. Vershynin. No-gaps delocalization for general random matrices. Geometric and Func- tional Analysis, 26:1716 1776, 2016.
Jun Shao, Yazhen Wang, Xinwei Deng, and Sijian Wang- Sparse linear discriminant analysis with high dimensional data. Annals of Stalistics, 39:1241 12622111265, 2011. doi: 10.1214/10-A05770
Houssem Sifaou, Abla Kammoun, and Mohamed-Slimn Alouini. High-dimensional linear discriminant analysis classifier for spiked covariance model. Journal of Machine Learning Research, 21(112):1-24 2020a.
Houssem Sifaou, Abla Kammoun, and Mohamed-Slim Alouini. High-dimensional quadratic discriminant analysis under spiked covariance model. IEEE Access, 8:117313 117323, 2020b.
Tgor Silin and Jianqing Fan. Canonical thresholding for nonsparse high-dimensional linear regression. Annals of Statistics, 50(1):460 486, 2022.
R. Tibshirani, Trevor Hastic, B. Narasimhan, and Gilbert Chu. Diagnosis of multiple cancer types by shrunken centroids of gene expression. PNAS, 99:6567 6572, 01 2002.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statislical Society, Series B, 58:267-288 01 1996. doi: 251111111117--161199999.6
Ryan J. Tibshirani and Jonathan Taylor. The solution path of the generalized lasso. Annals of Statistics, 39:1335 1371, 2011. doi: 10.2214111-OSS1
Roman Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cam- bridge University Press, 09 2018. ISBN 9781108415194. doi: 10.1017/97811002331599
Vincent Vu and Jing Lei. Minimax rates of estimation for sparse pea in high dimensions. In Proceedings of the Fifleenth internalional conference on artificial inlelligence and statistics, pp. 1276 1286, 2012.
Weichen Wang and Jianqing Fan. Asymptotics of empirical eigenstructure for high dimensional spiked covariance. Annals of Statistics, 45:1342-1377 1374, 06 2017. doi: 10.1214/16-A0S4487
Daniela Witten and Robert Tibshirani. Penalized classification using Fisher's linear discriminant. Journal of the Royal Statistioal Society, Series B, 73:753-772, 11 2011. doi: 10.1111111467-998820011 0tt.-99888.21100775800 2011.0078..x
J. Yao, S. Zheng, and Z. Bai. Large Sample Covariance Matrices and High-Dimensiond Data Analysis. Cambridge University Press, 2015.
Peng Zhao and B. Yu. On model selection consistency of lasso. Journal of Machine Learning Research, 7: 2541 2563, 12 2006.
W. Zhao, R. Chellappa, and A. Krishnaswamy. Subspace linear discriminant analysis for face recognition. In Third IEEE conference automatic face and gesture recognition, pp. 1998.
23

---

Published in Transactions on Machine Learning Research (January
Mu Zhu. Discriminant analysis with common principal components. Biometrika, 93(4):1018-1022 2006.
Yinchu Zhu and Jelena Bradic. Significance testing non-sparse high-dimensional linear models. Electronic Journal of Statistics, 12(2))3312 3364, 2018, doi: 10. 10.1214/18
EJS1443. 10.114/18-E.S11443
Verena Zuber and Korbinian Strimmer. Gene ranking and biomarker discovery under correlation. Bioinfor- matics, 25:2700-211 2707, 08 2009. doi:
bioinformatics/btppoe.
7 Supplementary Numerical Results
7.1 Supplementary Results to Model 1 to Model 3
We provide here several tables that supplement the simulation results in Section I Table 9 reports the magnitudes for the strong and weak signals in the discriminant dircction and whitened direction Ç for Model 1. Table 10 and 1 reports the average number of strong signals and weak signals captured by Ida pea, DSDA, AdaLDA and SLDA under Model 1. From Table I we that Ida pea includes only few features corresponding to the weak signals and thus selects only a small number of features compared to the other classifiers. We also report in Table 12 the summary statistics for the condition mumbers of the covariance matrix Σ in Model 3; these statisties indicate that the performance of AdaLDA for Model 3 can sub-optimal as AdaLDA estimates 3 by solving a linear programming problem which numerically unstable when Σ have large condition numbers.
Table 9: Strong and faint signals for equal correlation setting (model 1). See the discussion Remark
P
0.50
0.60
0.70
0.80
0.90
strong signal /
1.98
2.47
3.29
4.94
9.88
strong signal Ç
1.40
1.56
1.80
2.21
3.12
faint signal /
-0.02
-0.03
-0.04
-0.06
-0.12
faint signal Ç
-0.02
-0.02
-0.02
-0.03
-0.04
Table 10: Average number of strong discriminative and whitened variable with standard deviations in paren- theses for equal correlation setting (model 1)
0.50
0.60
0.70
0.80
0.90
Ida pea
9.59 (0.82)
9.42 (0.85)
8.53 (1.41)
8.53 (1.41)
3.68 (0.98)
DSDA
10.00 (0.00)
9.99 (0.07)
9.99 (0.07)
9.99 (0.07)
10.00 (0.00)
AdaLDA
9.95 (0.22)
9.93 (0.26)
9.95 (0.23)
10.00 (0.07)
10.00 (0.07)
SLDA
9.98 (0.22)
9.98 (0.28)
10.00 (0.00)
10.00 (0.00)
10.00 (0.00)
Table 11: Average number of weak discriminative and whitened variable with standard deviations in paren-
Cheses for equal correlation setting (model 1)
0.50
0.60
0.70
0.80
0.90
Ida pea
2.44 (4.22)
1.90 (3.73)
0.98 (3.35)
0.98 (3.35)
0.00 (0.00)
DSDA
86.20 (31.06)
96.42 (32.57)
107.16 (30.53)
107.06 (29.96)
86.48 (8.66)
ADaLDA
36.59 (5.94)
35.91 (5.22)
36.78 (6.36)
37.86 (5.88)
38.65 (5.02)
SLDA
718.86 (226.57)
778.18 (95.53)
789.96 (0.21)
789.96 (0.21)
789.96 (0.21)
24

---

Published in Transactions on Machine Learning Research (January /202)
Table 12: Summary statistics of condition numbers of generated Σ from uniform, normal and student distributions
Statistics
U(-1,1)
N(
ts
Mean
415.3
688.0
1059.2
SD
102.4
195.2
400.0
Median
391.5
658.5
966.3
IQR
108.0
205.9
314.1
Max
849.9
1622.1
3708.4
Min
266.3
396.9
536.0
25

---

Published in Transactions on Machine Learning Research (January /202)
7.2 Sensitive Analysis
This subsection aims to explore the sensitivity of Ida pea to the number of spikes d. As a baseline, we list the performance of Oracle and DSDA. Consider a covariance matrix Σ = (ox)px with parameters P1, > 0 and P1 + P2 < 2:
= uupuiu + Pupuzu + (2--p1 - -P2)(Ip - -n u2u?) where U = vp 1p and U2 1/2 (e) (e)- e2).
(7.1)
This covariance structure is slight modification of Model (the equal correlation model) to incorporate two spikes corresponding to the eigenvalues PiP and P2p. The class means set to be μ4 Op and C1 + 20, so that the mean difference #2 P1, discriminant direction 3 and whitened direction Ç are all sparse with exactly three nonzero entries, located at the first three coordinates. This contrasts with the dense structures exhibited Model and Model 3.
For sensitivity analysis, our method variants are labeled as lda o pea(d) for 1,2,3, with being the correctly specified case. Training and test samples are generated using from Eq. 77.1. following Monte Carlo setups and CV procedure detailed in Section 4.1 Table 1 shows that correct specifying overestimating results in classification errors comparable to DSDA whereas underestimating d leads to notable degradation in performance. Since overestimation does not harm accuracy, we recommend, in practice, selecting as the smallest integer such that the cumulative eigenvalue ratio A))tr(() 2 0.9, ensuring 90% of the total variance is retained in the selected subspace.
DSDA achieves similar accuracy but tends to over-select non-informative features as shown in Table with average 49.46 selected features well above the true sparsity level of 3. In comparison, Idaopca maintains comparable accuracy across all simulation scttings. Our method also yiclds more parsimonious and offers greater compulational efficiency especially under the spike covariance structure.
Table 13: Mis-classification rate (%) with standard deviations (%) parenthescs for Eq. 7.1 based on 200 independent Monte Carlo replicates.
(21.P2)
(1.5, 0.3)
Oracle
0.30 (0.36)
Ida o pca(d 1)
24.19 (23.09)
Ida o pea(d 2)
0.55 (0.60)
Ida o pca(d 3)
0.55 (0.60)
DSDA
0.50 (0.49)
7.3 Elliptical distributions
A random vector E RP is said to have an elliptical distribution with mean j and covariance matrix Σ if the probability density function for X of the form
= -/2 - TI- p)), for all x E R..
(7.2)
Table 14: Average number of nonzero coefficients with standard deviations in parentheses for for Eq. based on 200 independent Monte Carlo replicates.
7.1
(p1.P2)
(1.5, 0.3)
Ida o pca(d 1)
1.49 (1.90)
Ida o pca(d -
3.82 (4.28)
Ida o pca(d 3)
3.71 (3.92)
DSDA
49.35 (32.97)
26

---

Published in Transactions on Machine Learning Research (January /2026)
Here Cy is normalization constant, |Σ| is the determinant of E and h is monotone function on oo). The class of elliptical distributions is analytically tractable and many results that hold for multivariate normal distributions can be extended to general elliptical distributions. In particular, Fang & Anderson (1992) showed that Fisher's LDA rule in Eq. (1.1) is also Bayes optimal whenever the feature vectors are sampled from a mixture two elliptical distributions with known covariance matrix Σ and known class conditional means μ1 and μ2- Shao et al. (2011) and Cai Liu (2011) leveraged this fact to show that, under certain conditions on the sparsity of either μ4 - μ22 and Σ or the sparsity of Eー -2=), both SLDA and AdaLDA classifiers also achieve the Bayes error rate for elliptical distributions with unknown Σ. Ma and P2. We now discuss how the theoretical properties for Ida pca can be extended to elliptical distributions, provided that they have sub-Gaussian tails as define below.
Definition Let V: [0, oo) + [0, 00) be non-decreasing and non-zero convex function with (0) 0. Let Z be mcan 0 random variable. The Birnbaum- Orlbauu V-norm of Z is defined as
|| inf { 2 EV EV
7 s 2
(7.3)
Similarly, if Z is a mean 0 random vector taking values in IP then its v-norm is defined as
7
sup WEBP, w
1 Z
(7.4)
Let v (x) exp([=:) and (2(2) = exp(=²). A mean 0 random vector Z is said to be sub-exponential if I v and is said to be sub-Gaussian if || [v, co. If Z E I2 is mean 0 sub- nub-Gaussian random vector then w Z is a sub-Ganssian random variable for all w E RR. Furthermore, if Z is a sub-Gaussian random variable then there exists a universal constant K such that for all & > 0, we have
I2 > t) s 2exp(=t
For more on sub-Gaussian random vectors, see Vershynin 2018 Section 2.5, Scction 3.4).
Remark 7. If Z E IRP is mean sub- Gaussian random vector with covariance matrix then || > w Eu for all w E RR. In this paper we shall assume that a converse incquality also holds, namely that there exist a constant 6 > 0 such that,
w T Ew 2 cillw T Z for all w E IR.
(7.5)
We note that the constant C1 in Eq. 7.5
can depend on Z but docs not depend on the choice of w E R._ If Z is multivariate normal then Eq. 7.55 a always hold. If is a zero mean sub- Gaussian random vector but not a multivariate normal then Eq. terms of its variance. This then allow us to obtaín better estimate for Σ Σ-> Σ in spectral norm, especially in the setting where the spiked eigenvalues could diverge with p. See for example Theorem 4.7.1 in Vershynin (2018) and Theorem 9 in Koltchinski & Lounici 2207
77.5 allows us bound the Orlicz norm of w TT Z for any w E RRP
We can now reformulate the classification problem in the earlier part of this paper to the case of elliptical distributions as follows. Let (EII...., and fe21...., E2ny} be independently and identically distributed mcan 0 random vectors with probability density functions of the form in Eq. (72)) and supposc that the training sample is given by Xij = #4 + tij for E (1,2) and E n1:).
Given these training samples { }. let X be the (n1 + 122) x p matrix whose rows are the {Xar). Then Theorem 1 in particular Eq. (3.4). also holds for those X as long as the {ay} satisfy Eq. (7.5. This then implies that Theorem 2 in particular Eq. (35), also holds when the X are sub-Gaussian random vectors. The resulting bound for - allows us to recover St by thresholding č. and hence Ridaooea + Rp. In summary, Ida o pea is asymptotically Bayes optimal whenever the feature vectors {X,}} are sampled from a mixture of elliptical distributions with sub-Gaussian tails.
27

---

Published in Transactions on Machine Learning Research (January /2026
8 Proofs of Stated Results
This section contains the proofs of Theorem | general assumption that the feature vectors {Xiy} are sub-Gaussian random vectors; see also the discussion in Section
through Theorem 3 We will present the proofs under the more
7.3
8.1
Preliminary results
We start by listing some elementary but useful facts about the 2 - oo norm and its relationships with other matrix norms. Recall that -| denote the spectral norm if its argument is a matrix and denote the la norm if its argument is vector.
Proposition 1. Let A Reixes and B E RR be arbitrary real-valued matrices. Let IR be arbilrary vector. For given i E [pi], let e; denole the ith elemenlary basis vector in Re Then
||A
max ||ATe:l) iE(m)
ezll;
(8.1)
||AI|1 $ 11-A112+000 ||>||:
(8.2)
|| S ||All2-ss0 x ||]||
(8.3)
Eq. of A. Eq. 8.2 provides bound for ||Axi]oo that is lighter than the naive bound || S || Axxl| < ||||]] ||]]]
881 states that the two-to-infinity norm of matrix oquivalent to the mazimum & norm of the rows
Throughout the section, will make of Bernstein's inequality, For completeness, state the result below without proof.
Proposition 2. Let X1, X2,..., Xn be independent, mean-zero, sub-exponential random variables. Then for all t > 0, there erists a constant e > such that:
x 21 I 2exp (-c- min
12 Σ 1= 1
max; | )
P( 1=1
-ii
(8.4)
where || ll,, denotes the sub-ezuonneeee sub-ezponential Orlica norm.
Let E; = X; - p4. Recall that M) = (>x + 0 1/2 and 0² p- (tr(()r - tr(A)) are the eigenvalues of W = E 1/2 1/2 The following lemma provides several concentration inequalities for & E, îT - 77 and d² 6². Lemma Assume that the random variables {X,} satisfy Eq. (7.5) and the covariance matrix Σ satisfies A Assumplion and Assumplion 4 Then the sollowing bounds hold simullaneously with probabilily al least. D x-2 (where En is defined in Eq. 2.4h)
1ΣS Σ|| O(p n In
(8.5)
1|E.11 O(p - Inp),
(8.6)
Σ Σ OOpv n Inp).
(8.7)
fmu mk O( 1n Inp), for all k E [d]
(8.8)
2 0 O0-nn1nn1p
(8.9)
Eq. 8.5 is given in Lounici 2014) and Koltchinskii & Lounici (2017) while Eq. 8.6 follows from an
application of Bernstein inequality; see Section 2 and Section 3 Vershynin (2018) 1. Eq. 8.7 follows from
Eq. (.55 and Eq. (8.6) together with the observation that E E n (mi51 + 2 Finally, Eq. 88.
and Eq. 88.h follow from Eq. (8.7 and Weyl's inequality.
8.2 Proof of Theorem
Recall that 24 and á denote the px matrices whose columns are the orthonormal eigenvectors corresponding to the d largest eigenvalues of E and the pooled sample covariance matrix E, respectively. Now let U and
28

---

Published in Transactions on Machine Learning Research (January /002)
Ú1 be the p (p- d) matrices whose orthonormal columns the remaining cigenvectors of E and Σ, respectively, i.e, Ip uu = U.Ms and L, -uũ = =ū
Let =\t be the d x d orthogonal matrix that minimizes
min ||W - U úl1F W
among all orthogonal matrices. Let En - Then by Theorem 3.7 in Cape et al. 2019 we have
| - 4- s 2(Aa 02) - (2, )
+ 2(A + - (E, (2(22) 1112-200 x || sis 20)
(8.10)
+ 2(A + 12) )I(4) x sin ((2,24))
+ || sin ((2,22))2 ||
Now recall the matrix Eo from Eq.
2.4
We then have
Σ
- E
(I1
(1)
11)
n22
(32 12) 22)
8m
F..
E
E2
Using the same argument that for the proof of Theorem 1.1 in Cape et al. 2019 have with probability at least 1 - that
1/2
(A1 + o²) Inp
, E 112-000 S Cd max Eji iElp)
x
(8.11)
(d1 + 02) Inp
u E. S COO/
(8.12)
n
where we have used the assumption that A ((p) for all k E [d] so that r(2) the effective rank of Σ is bounded. Here and in the subsequent derivations we will, for simplicily of presentation, use to denote a finite and universal constant that can change from line to line.
Therefore complete the proof of Thcorem ] suffices to show that, for j {1,2), the terms | Ui E, E. 1|2-100 and (NV2TE,A 22 of the same or smaller order than those in Eqs. 8.11
and 8.12 respectively.
We now bound
11a-00.. From Assumption 4 and Proposition we have
Iuuul S <Cd, and E E,U 1120000 I [Vd
(8.13)
Furthermore, we also have
|E S Vd max VE(E
eE
(
(p) 13}l = Vd
max (e(r)e(dd
-11 x [[X m1) uu;l-
Since X, is sub- Saab-Gnsssinn, by the properties Orlicz norms we have
I1((7 11) T Te! [(X n1) u,11l < - u1) Te x || m1) T Tullll
Eq. 7.5
implies that there exists a constant C > 0 such that for any E (x] and j [d].
||x; m1) of s sc
In C
(max Ens
1/2
(8.14)
nn
/ni
(Elp)
+
|[x, - -11 T (ullo, Var(u] (X1 -1) s C n -Eu, S sc
(8.15)
n1
29

---

Published in Transactions on Machine Learning Research (January /202)
Now fix arbitrary pair (i,j) with i e pp] and j [d]. Then by Assumption 4 Eqs. 88.4 properties of sub-exponential random variables, we have
and 8.15 and
e (P) ا,ررl C(11 + o²)
d
P
Furthermore, by Bernstein inequality with probability at least 1 - -(ppp))
/Vershynin 2018
Section 2.8), there exists constant C > 0 such that
(p) ,2)1 K
E u,)l + n1
(
Clup
max Eai
1/2 VAi + o2.
(8.16)
[E(p]
Now recall Eq. (8.1)). Then by Eq. 8.16 together with union bound over all i [p] and j [d] we have, with probability at least 1 Opp ²),
E C S Cd max (Ere u,yl = C d3/2 Vplng
(8.17)
(.)
n
where we had used Assumption namely n, and A x P. when simplifying the above expression. An almost. identical argument also yields
14 Eguu SC" //pppppp Vplup
n
with probability at least 1 -pp-0ppp))
We therefore have
2 E, 1120000
C
(d(maxX<>) Ex) 1/2
12 x
(A1 + Inp
1 In Vp
Ad 1 o²
As +
o²
7
(8.18)
Inp
CdB/2
S
np
We next consider |14 U[ E, 12
For j d, we have
[I(X; m1) uull S
S Cn 1/2 o.
Then following the same argument as that used for showing Eq.
8.17
we have
1/2
ln S
ln
11 112 S Vpdo max E.s)
CCo/pd
(8.19)
(E(P)
n1
n
with probability at least 1 Opp-2), and similarly for |0 E U 12-000
Next have, by the Davis-Kahan theorem and Lemma I that
||
Inp
|| sin 0(Ū,2))) < Aa o²
S
(8.20)
n
with probability at least 1 - Opppp 2). From Eqs. 8.12 1.19
and
8.20
together with a similar argument
as that for showing Eq. 8.18 we have
2 u E
112400 sin O(Ú,4)))
Inp n/P
& Cod x
(8.21)
(Aa+
with probability at least 1 - O(p Eq. 8.20 together with Assumption also imply
sin 1(1 x | Cluy Vd Cvdlnp
(8.22)
n/p
30

---

Published in Transactions on Machine Learning Research (January /2026
with probability at least 1 00ppp- Next note that
|u 12- = | u1 11 S 1 =
We therefore have
inni _u_
2000 x
sin ((2,21)
Co²
Inp
Co VIp Vnp
(8.23)
Aa +
Ad + o²
Subtituting the bounds in Eqs. 8.18
8.21
8.22 and 8.23
into Eq.
8.10
we obtain
| - us 2-000
S
c(
d/2
Vinp
In p
+
o² V'pp
sc
d Inp
Vnp
nyP
Vup
with probability at least 2- This completes the proof of Theorem
8.3 Proof of Theorem 2
First recall that ? = WX2x X1) and Ç - W(22 - μ1) where the whitening matrix W and its estimate W
are given by Eq. 3.2 and Eq. (3.3), respectively. We now consider the decomposition
s-c= w (X2 X1) (m2 - -1)] (W () X1) (u2 m1) + (W WW) m1):
We will now bound each of the term in the right hand side the above display. We start with the term in (A). Let 8 (X2- X1)= (m2 -u and let { = W8. We then have E[g] and Var[E] c 'I, for some finite constant Since ā satisfies Eq. (7.5,, { also satisfies Eq. (7.5)- Hence, by Bernstein inequality for sub-Gaussian random vectors, there exists constant C > 0 such that with probability at least 1 - 0(p
[ X1) - (#2 - (1) = 11 S sc /c/"/
mm
(8.24)
We now bound the terms in (B) and (C). Let D and D be diagonal matrices where
D = (Å + I2)
1/2 1/2 D= oZI)) 1/2
We start by decomposing w W as
W- W --ũ UDU T + 2 uu + 5- --nnT - -UU
(8.25)
111
Now consider the term (o - - LuT)8 obtained by combining the expressions in (B) and (IT). The covariance matrix for 6 is n n-E and hence (Ip - UU) 8 satisfies Eq. 7.5 with covariance matrix n (I, UU ). Therefore, by Bernstein inequality, there exists a constant such that with probability at least 1 - O(p ?),
Inp I|ZZ, -uu -uu )8 CCO
Furthermore, from Lemma 1 we have with probability al least. 1- O(p 2) that
102 - o² 80(0 + +
lup
10 18-1 10-100-1-
S SC
(8.26)
Combining the above bounds, we obtain
Clnp
((6-1-0-11 uu )81 10T v
(8.27)
n
31

---

Published in Transactions on Machine Learning Research (January /2026
with probability at least 1 1 Opp
We next consider the term ( "Ip - uu uu (2z v1). Recall that if 4 has bounded coherence
in Assumption then |[ - 1444111a I (1 +u)vv where Cu is finite constant. We therefore have, by
Lemma that
|1(6 - o "I - uul ( )2 u1) 1 0 ([ - Uu (u2 Mu2 pi)
|[ o "I(Tpp - -ut) 21
160 -1 - x o(1 + Culva 115
lnp
Cy
1|01 n
with probability at least 2).
We now focus our efforts on terms involving ÜDůT TDUT Let E be the minimizer of I| u uūul among all d orthogonal matrices W. We then have
ûBũ -DD Uyyu + LuTúbbū UDU
(ύ - (1) (3))
-uu ú =) =)
a DŪ uu AuDi UDU
Part 2
Part
We now note a few elementary but useful algebraic facts frequently used in the subsequent derivations.
Fact
NU - 311 < sin 201,21)1
(8.28)
jurTu || sin 0(4,20)) S |1a0T -uu 2|| sin @(Ū,2))) 2111
= 112T UI = u 4I1 = sin 0(2,20))
(8.29)
(8.30)
[ urūuD- DU' Tul & 114 (2 - (x-E x |||]
(8.31)
where (Hl,) is d d matrix with entries
Hig
V +0)+
+ Sr + 02)
Eq. 8.28 from Lemma 6.7 in Cape et al.
2019 while Eqs. 8.29 and 8.30 are standard results for
the sin-G O distance (see for example Lemma observation
Cai & Zhang (2018)) Finally Eq. (8.3)) follows from the
(u (uÚD-DuTu, = (U (UT D;:))
(Ai+ - (2g + +2)
= =(UT(Ois) (2, + ( ((), +02)1/2 + (A. + 02)1/2)
-AUTT - 41 (UA),,HH. (uT(T - H,
We thus have 4TuD D - H. Therefore, by Schur incquality for Hadamard product
(see e.g., Theorem 5.5.1 of Horn & Johnson
1991
have
I1(UT (I - (x- o H| & I1 ( EYHI x |
We next state technical lemma for bounding several terms that appears frequently in our analysis.
32

---

Published in Transactions on Machine Learning Research (January
Lemma 2. Suppose that Assumption 1 through Assumption 4 satisfied. Then with probobility at least
O(p-2), 2), the following bounds hold simallaneously
[ = p 3/2),
(8.32)
II - -E O n Inp),
(8.33)
lIeT ED2U1 = O n (n-1))
(8.34)
U 1uúD- - -D = Onn 1/2p 1/2 Imp),
(8.35)
NiăT - uu O ((vnnIlnp))
(8.36)
Eqs. (.33) and 8.34 follows from the sub-muliee sub-multiplicativity of the spectral norm together with bounds for
[ Σ from Lounici (2014) and Koltchinskii & Lounici (2017). Eq. (8.36) is a consequence of Eq. 88.30
and the Davis- Kahan
Kahan theorem. Eq. 18.32 follows from Weyl's 's inequality and the bound for -|- Finally,
Eqs.
(8.31). (8.22] and [8.34] together imply Eq. 8.35
With the above preparations in place,
resume our proof of Theorem 2
We first have
1(u - uu uyy-u 0l S (I0
1 vN + 1 |
=||| 8|l
4 C(1
d3 In p
vdlmp)
|iBi 811
n/p
d3 In p
vālnp) p)
((12] + +ID DDDDI 81
n/pp
(8.37)
d3 Inp
vdlnp)
Inp np
|1 ((iTTo) 61l
+ "
n//P
Vdi +
Cylp
|1 [2TT8||
PVa
with probability at least 1 1-0 For the above inequality, we have used Lemma to bound _
- DII and
used Theorem to bound 112-- - 2- Finally we used Eq. 8.281. Eq. 88.30 and Eq. (836 to bound lú El|.
Next let T =Ip Vu, Then
-DDbb I I UD -DU -DUTUÇÇT 811 + |ID TT8|
|111225000 (1121 uTD - DU Ôll |D|| x 1 T
(8.38)
(vaamp + Cvalup) x plnp Pvn
Clap
Pvn
nvp
with probability at least - 0pp-2). In the above derivations. we bound (TT - Du using Eq. 8835
and bound |/2 T |44 (Ip ---l -niT) ]l using Eqs. 8.29 8.30 and 8.36 The bound for |D]] and ||2 ||
follows from Assumption and Assumption respectively.
Combining Eqs.
8.37 and 8.38 we obtain
N@ÜDūT uDDT)oll 0 Inp
(mm)) mppp
(8.39)
with probability at least 1 O(p-2).
33

---

Published in Transactions on Machine Learning Research (January /2026
Table 15: Asymptotic Order of Each Term:
Expression
E
Corresponding Terms
||Vllb
V
In
p
n.a.
(A)
-ÚDũ -UDU -uDu T
In p nvP
p
|||||
(B)-(1) and (C)-(I)
np
1160-1 - --- -UM -ulT)vlls 110 100- - -uuT)ala
In
V
lu
||
(B)-(IT) and (C)-(TI)
In
In
||c||
(B)-(7IT) and (C)-(IIT)
Using similar arguments as that for Eqs. 8.37 and
8.38
we also have
|( - uE) - uu =))x 81100 -Iuuu To
lini uuu))ll S
d3 Inp
Vdlnp
plup
dlup np
plup
<c(1
V
+c
V
(8.40)
np
n/pp
0 (np)
with probability at least 1 - Opp 3), where once again T Ip ÜÚT Note that in the above derivations have used Lemma ! to show that 8-1 1 is bounded away from 0 by some constant not depending p and n.
Finally we can also replace with g in the derivations of Eqs. 8.37 1,
8.38
and 8.40 to obtain
I(ĀDū - -UDU ) =
In (
(8.41)
npp
||0 | -uu ) =0 Inp [KClla
(8.42)
n
simultaneously, with probability at least - O) A summary of the bounds for the terms (A), (B)- through (B)-(IIT), and (C)-(T) through (C))(III), are provided in Table 15 Combining the terms in this table we obtain the bound for I| - <lloo given in Eq. (3.5) (note that both ||||| and [|C]] bounded, see Assumption 2 This concludes the proof of Theorem
Finally, for ease of reference, we state two collaries summarizing the main derivations in the proof of Theorem These corollarics will be used the proof of Thcorem below.
Corollary 1. Suppose that Assumption 1 through Assumption 4 are satisfied. Let v be either a fixed vector in RR or a p-variate sub- Gaussian random vector with E(v] 0. We then have
(n max; 5t) (n 1/2 1-12///
if v mub-Gaussian vector
=
if v is constant vector
with probabilily al least 1 - O)) where s? is the variance of the tth element of Corollary 2. Suppose that Assumption through Assumption are satisfied. Let be p-variate sub- I for some finite O. We then have
Gaussian random vector with Var[v] = c n
=
with probability at least 1 O(p-2).
34

---

Published in Transactions on Machine Learning Research (January /2026
8.4 Proof of Theorem 3
For simplicity of notation we will write S instead of Se since Theorem 3 Bonly depends on the whitened veclor Ç. Now recall the definition of S as
$ {j 18g1 tn where tn = [mp for some constant < a
< 2.2 We now show $ = S asymptotically almost surely.
From Theorem there exists choice of C such that if Bn CV n Imp then
p ( &1 B.1) 3)
j=1
Now suppose 8 n S * 0 where (-) denote set complement. Then there exists such that S 0 and GG)= ba and for this j we have KG را > tn > Bn, provided that n sufficiently large. We thus have
PPSS x S P(U &1>> 11) 0p
(8.43)
Similarly, S So + 0 then there cxist a j such that 1cyl > C. and (Ś91 S tn- Recall that Co > 0 is the constant appearing in Assumption 2 in particular, Co does not depend on n and p. By the reverse triangle inequality, Ky; - رن Co t, > Bn for sufficiently large and hence
P(Sn Sc + 0) SP =SP(UUIS
- B..) 8..)= = =(p
(8.44)
1-1
Combining Eq. 8.43 and Eq. (8.44) yields IPSS + $) O and hence, as p -> 00, by the Borel-Cantelli lemma we have S asymptotically almost surely.
We now show that the error rate for lda pea converges to the Bayes error rate Ry asymptotically almost surely. From the description of lda pea in Eq. 2.9 it is sufficient to show that
C3 (WZ
81+ 82
)
<5 W[ #it P> 0.
(8.45)
P
In
(8.46)
72
1- क
Note that the convergence in Eq. 8.45) is with respect to a random testing sample TT, Σ) + (1 m1),(,, 2) together with the training data, while the convergence Eququatioo 8.46 is with respect to the training data only. As Eq. equation 8.46 follows directly from the strong law of large numbers, we thus focus our efforts showing Eq. (8.45).
23 First, suppose Š S and let h(Z)= = 23 (i - [W = Then
[M(Z]
so
[i(Z x1t881) W(z - m1ta2) ( )lll
+ so
[W(z = )} sll ul5 0l + so [i x1+82) w(Z
(8.47)
) 11 Clloo- uu 2
The bounds for o/ and [| &o given in Assumption 2 and Theorem respectively. It thus suffices to bound
(D) = || - x1t&1 - W(z - /1 )
(E) = [w / )s
Write Z = M4 + E1/2e, where μ4 = M4 if Z is sampled from class and μs = P2 otherwise. The term together
11OV - w can be analyzed using the decomposition for W W given in Eq. 8.25
35

---

Published in Transactions on Machine Learning Research Sanuaryy20200
with similar arguments to that for deriving Eqs. 88.7 probability at least. - O(n that
8..7
88.38
and 8.40
In particular we have, with
Inninp /
I[OW- ( O(V
(8.48)
n
Next, using Corollary 1
Corollary and Eq. 8.48 we obtain, with probability at least 1 OOn
(D) s [1(0 - W)E
ssllo 1100 W)
_ - Waill =110x1; - WM2lpa
s <clv Innlup
Inp
= Inp
In Inninp
n O(V
n
n
Thirdly, we have
(B) = - )l + [w -z) = 1151loo + +1leesslllo =: O(Z)
(8.49)
where [e=]s is a mean sub-Gaussian vector in & and Var[[ea]s] Izg. Therefore, by Assumption and properties of sub- Gaussian random vectors, the term (Z)) is bounded in probability.
Combining the above bounds, we conclude that with probability at least 1 - OOn 2), S S and
Di 2 x1 - c3
(ss &Coo (V
111100 + =d(Z)) Inp Innnnp
Hence, for lnp o(n), we have
Di -
- (3 - 1 ]s 0
in probability. This completes the proof of Theorem 3
8.5 Theoretical results for Section 5.1
We now present theoretical results for multi-clas. Ida pca. We first assume that the whitened directions ) the whitened indices {St}, and the covariance satisfy the following two conditions, which are natural
generalizations of Assumption 2 and Assumption 3 to the multi-class setting.
Assumption 6. Let |S; = > 0 for each i = 2, K. Recall that S is the set of indices j for which Siz # 0.
Let Co > 0, M > 0 and Cę > 0 be constants not depending on such that max, S M and
min min Giyl 2 Co:
max KE[K]
| 1/2 1/11 S Cc.
K JES,
Assumption 7. Let 0 be fixed and that, for sufficiently large p, 11. nk and p satisfy
T4 n
Inp onn). (1) for i j. i.ij [K] and
Furthermore, for sufficiently large the spiked eigenvalues Ai,..., da satisfy
λ O(p), for all ls E [d].
Given the above conditions, the next result. extends Theorem (and has an identical proof) to bound kÇ - -(0) ooo for i 2 2.
Thcorem 4. Under Assumption 1 Assumption Assumption C and Assumption 7
there crists constant
C > such that with probability at least Opp
max |1 - -((0) 11 S scv Imp
(8.50)
$=
36

---

Published in Transactions on Machine Learning Research (January
We now consider hard thresholding estimate (0) for recovering ( For given i 2 2, let
"- 2 1(1)) 1.),
E (p]
(8.51)
where i, (Iup/n) for some constant 0 < a < z.4 Here assume, for simplicity, that a takes the same value for all classes. Given ((t) define the associated active set S; = {j: (00) * O}. The next result extends Theorem 3 (and has an identical proof) to show that Ida o pea is also asymplotically Bayes-oplimal in the multi-class setting. However, we note that (to the best of our knowledge), there no closed-fore formm explicit expression for the Bayes error Ry when classifying data from a mixture K 2 multiviarate normals.
Theorem 5. Suppose that ~ Σ 1 TuN Σ) where nı 2 and Σ R, 1. Suppose Assumption
Assumption 4 Assumption 0 and Assumption are satisfied. We then have
max P(S. + S₂) Owp 3). 1=2, K
(8.52)
Furthermore we also have Ridaopca- RØ + almost surely as n, - 00.
8.6 Theoretical results for Section 5.2
We now present thooretical results for qda o pca. which depend on the following variant of Assumption 2 through Assumption | for heterogeneous covariance matrices.
Assumption 8. Let |Ai| = > 0. Recall that A; is the set of indices for which Gij + 0, = 1,2 where
Si 2, 1/2 p4. Let Co > 0, A and Ce > constants not depending on p such that max a1, a2 s M
and
min min (say 2 Co: (E(1.2) jEA;
max max (E(1,2) (6(12)
|Σ 1/2 u S CCC.
Assumption 9. Let o; be fixed and that, for sufficiently large p. nı, n2 and satisfy
f1 = 0(1,, Inp- = o(n).
Furthermore, for sufficiently large the spiked cigenvalues Aa....., Xid, satisfy
Ayk = Op),, for all k E [d.], i = 1,2.
Assumplion 10 (Bounded Coherence). There is a constant Cu 2 1 independent of n and such that
1411 S Cuvd, vD for i = 1,2.
Assumption S guarantees that the noncentrality parameters for these xi are strictly positive and finite, that the Bayes error rate Rr is stricly bounded away from 0 and min(TT, 41} (which corresponds to random guessing). Note that if Σ + Eg then there is no simple closed-form expression for Rp as it depends on the tail behavior (2003) Section 6.10) for more details. Assumption g allows the spiked eigenvalues for each Σ, grow linearly with the dimension p. in contrast the bounded eigenvalues assumption frequently encountered in the literature Li & Shao 2015 Cai & Zhang 2021
linear combination of independent, noncentral xi random variables; see Anderson
We then have the following extensions of Theorem and Theorem
Theorem 6. Under Assumption
through 10 there exists constant C' 0 such that
Inp n
max 11.. Sillo s e
(8.53)
with probability at least 1 -pp-2).
37

---

Published in Transactions on Machine Learning Research (January /2026
Construct &e as in Eq. 3.6 for {1,2} and let A be the indices for the non-zero coordinates of &t. Let Z ~ TING(11 En) + (1 (U, E2). Then
max P(A # A.) 3)
(8.54)
Furthermore, we also have Rydaopea R + in probability, as n, p- - 0c.
Proof. For conciseness, we omit the derivations of Eq. (8.53) and Eq. 8.54 as they follow the same argument
as that in the proof of Theorem and Eq. (3.)
In order to show that Q(Z X., Wi, Ao) is consistent estimate for Q(Z Mi, Wr, Ao), it suffices to show that the following quantilies
(Z) = 1 Z - (-:1))) Wi (Z - -1)1 M (z - -x)] M4 - -Xi]
(8.55)
(z(Z): = [[aa [2) [M2 ((0(Z-==))))) [2 - -xx) [22 -X2]]
(8.56)
both converge to 0 20.
Assume without loss of generality that Z is test sample from class 1, i.e. Z ~ NNp(11 E1). Then
(i) S a0 aol||Wi (Z m1) - VZ -11)
(8.57)
+ +2o0| W2 - [1)llall IIMi p1) 1Z X1)lAo
where 40 = Aol & a1 - a2 with a, = |Agl.
Let (Z) = [[ 11) Wi (Z X1)lAa 21 Then by Corollaries 1 and and following the same
arguments that for Eq. (8.48 have
(Z) < 1(W Wi)mi 1000 + [1((W] - W/ A + |M, - Wiai
(8.58)
<c(v
lup + Inppup
Inp
with probability at least 1 (2), where cs ~ NOO,I))) Ip). An almost identical bound also holds for lZ( = [ (Z - 12) (Z 82)1Aa booo namely
h2 () < |lW - W2) + +I( - -W3) ex + 1M - Wopelloo
(8.59)
<c(V
Inp Σ 1/2 M1 1|2 +
Innlup
Inp
with probability at least 1- O(n ?).
Next, 11 [ (Z (4))] Ao a1 and Wi 11)] A.ll I1. be bounded by the quantity ((Z) defined in Eq. (8.9 which, by Assumption g is bounded probability.
Finally, by Proposition and Lemma in Jiang et al. (2018), we have jÂ-ii -x] + 0 in probability. Combining the above statements yield the proof of Theorem
38


Reviews:
[{"review_id": "x3kKMYqwC1", "reviewer": "Reviewer_fxez", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "audience_yn": "Yes", "audience_explanation": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "requested_changes": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "broader_impact_concerns": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "additional_comments": "Please disregard this as my review was mistakenly made as an official comment, and my questions have been addressed.", "extra_scores": {}, "extra_text": {}}, {"review_id": "pvuk1O2M6w", "reviewer": "Reviewer_eQPd", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors study classficiation using the PCA + LDA approach, wherein first a PCA is used to decrease the dimension of the data and then LDA is applied for classification. They show theoretically that this procedure reaches the Bayes error asymptotically under mild assumptions, in particular only assuming that log(p)/n  goes to zero. They also extend their results beyong Gaussianity and two classes.", "claims_and_evidence_yn": "Yes", "claims_explanation": "I am not an expert on the underlying theory and literature, but the results and proofs, as well as the empirical results, appear solid. The whole paper including the proofs are well-written and the reasoning clearly explained.", "audience_yn": "Yes", "audience_explanation": "I think these are genuinely interesting results for a broad audience in statistics and machine learning. I was not even aware that PCA + LDA are used together and the theoretical results the authors present give a strong argument for this procedure.", "requested_changes": "This paper was a pleasure to read; well-written, well explained and with relevant results and an impressive literature overview. The proofs are well-written as well and easy to follow, with several smart writing choices that help the flow of reading, such as collecting often used algebraic facts in one place. \n\nMy only real \"concern\" is the length of the paper, specifically Section 5, which not only makes the paper longer but also disturbs the flow of the story. The way I understand the contribution, even without Section 5 the paper would be interesting enough for TMLR, so there is the question of whether this section is really needed. On the other hand, I see that that might not be enough for a separate paper, and it may also not be ideal to move it into the appendix.\n\nAside from that I only have minor comments/requests:\n- What exactly are the assumptions of d and s relativ to p? In the experimental section, Model 1 has no sparsity in the withened data, i.e. s=p. Is this a problem from a theoretical point of view or is this allowed? If so, would the method correctly identify p=s as the sample size increases? This may be spelled out more clearly in the main text.\n- This is a detail but it may be worth writing that Step 1 on page 5 is the same as PCA, depending on the audience this might not be immediately clear.\n- Page 11, when talking about Sifaou et al relying heavy on random matrix theory, maybe it could makes sense to very briefly spell out again what you do differently to allow for $p/n \\to \\infty$.\n- The Bernstein Inequality for random vectors is used several times in the proofs. It might be nice to quickly add it somewhere in the beginning of the proof section.\n\n\n- page 2 first paragraph, \"features selection\" --> \"feature selection\"", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "AD3sQlcSvI", "reviewer": "Reviewer_2mqw", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces and analyzes a classifier named lda$\\cdot$pca for high-dimensional data, where the number of features p can be much larger than the number of observations n. The core idea is to first use Principal Component Analysis (PCA) for dimension reduction and data whitening, followed by applying Linear Discriminant Analysis (LDA) in the reduced and whitened space.\n\nThe algorithm itself is natural and straightforward, while the main contribution of this paper lies in the theoretical part. The authors provide a rigorous theoretical analysis, demonstrating that under certain conditions, the lda$\\cdot$pca classifier is asymptotically Bayes-optimal. This means its misclassification error converges to the lowest possible error rate as the sample size increases. This is a significant result, as it is claimed to be the first of its kind for a classifier that uses PCA for dimension reduction prior to classification in this high-dimensional setting. \n\nWe also want to mention that the theoretical framework is notably built upon a spiked covariance model where the large \"spike\" eigenvalues can grow with the dimension p.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The claims made in the submission are, for the most part, supported by accurate, convincing, and clear evidence.\n\nClaim: The lda$\\cdot$pca classifier is asymptotically Bayes-optimal.\nEvidence: This claim is supported by Theorem 3, which states that under Assumptions 1-4, the misclassification rate of lda$\\cdot$pca converges to the Bayes error rate.\n\nTheorem 1 provides a crucial bound on the estimation error of the eigenvectors of the covariance matrix, and Theorem 2 bounds the estimation error of the whitened direction ζ. These intermediate results build the foundation for the main theorem. The authors also clearly distinguish their theoretical assumptions (e.g., diverging eigenvalues) from prior work, strengthening the novelty and relevance of their result. The evidence is convincing from my perspective.\n\nClaim: lda$\\cdot$pca is competitive with state-of-the-art methods in terms of accuracy.\nEvidence: This is supported by the numerical results in Tables 1, 3, 5, 7, and 8.\n\nClaim: lda$\\cdot$pca selects a smaller number of features.\nEvidence: This is supported by the model size comparisons in Tables 2, 4, 6, 7, and 8.", "audience_yn": "Yes", "audience_explanation": "Yes, from my perspective, the findings of this paper would be of interest to a significant portion of TMLR's audience. Here are the primary reasons:\n\n\nHigh-Dimensional Data: The \"p >> n\" problem is a central challenge in modern machine learning. Any novel method that provides a practical solution with strong theoretical backing is highly relevant.\n\nTheoretical Guarantees: The machine learning community, and TMLR's audience in particular, values not just algorithms that work well in practice, but also a clear understanding of why and when they work. The paper's main theoretical result—Bayes optimality under a realistic spiked covariance model—is a strong and appealing contribution.\n\nFeature Selection and Interpretability: The paper's finding that lda$\\cdot$pca often selects a much smaller feature set is a major selling point. In many applications, building a simpler, more interpretable model is as important as achieving high accuracy. This aspect of the work will appeal to researchers and practitioners focused on explainable AI  and model parsimony.", "requested_changes": "The paper is well-written and the contributions are significant. Following are some of my personal thoughts:\n\n1. The classifier has two main tuning parameters: d (the number of spikes) and s (the sparsity level of the whitened direction). For d, the paper suggests choosing it to account for 90% of the variability. For s, it is selected via cross-validation. While the paper mentions that \nd can be consistently estimated, a more detailed discussion or sensitivity analysis regarding the choice of the 90% threshold for \nd would be beneficial. How sensitive is the final performance to this choice? A brief experiment showing the impact of varying this threshold could add practical value.\n\n2. Computational Complexity: The paper does not include a discussion on the computational complexity of the lda$\\cdot$pca algorithm compared to the other methods. Adding a brief section or paragraph on this would be helpful for practitioners who need to decide which method to use on very large datasets. The PCA step is **generally efficient**, but a formal comparison would be a valuable addition. (There is no need to add computational complexity comparison in the review round, just my curiosity.)\n\n\n3. Performance on Dense Signal Models: In the simulation studies for Model 1 (equal correlation) and Model 3 (random correlation), the true whitened direction ζ is dense (all entries are non-zero). The paper notes that lda$\\cdot$pca works well by capturing the few \"strong\" signals and ignoring the many \"weak\" ones. While this is a strength, it would be interesting to see how it performs in a scenario with a truly sparse ζ that has only a few non-zero entries of similar magnitude. This would be a direct test of the feature selection mechanism in a less ambiguous setting and would complement the existing Model 2 (block diagonal), which is sparse in both the original and whitened domains.", "broader_impact_concerns": "No.", "additional_comments": "No.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)