A New Forward Discriminant Analysis Framework Based On Pillai's Trace and ULDA
Siyu Wang


Article Title]A New Forward Discriminant Analysis Framework Based On Pillai's Trace and ULDA










[1]Siyu Wangswang739@wisc.edu 







*[1]Department of Statistics, University of Wisconsin-Madison, 1300 University Ave, Madison, 53706, WI, USA









Linear discriminant analysis (LDA), a traditional classification tool, suffers from limitations such as sensitivity to noise and computational challenges when dealing with non-invertible within-class scatter matrices. Traditional stepwise LDA frameworks, which iteratively select the most informative features, often exacerbate these issues by relying heavily on Wilks' Λ, potentially causing premature stopping of the selection process. This paper introduces a novel forward discriminant analysis framework that integrates Pillai's trace with Uncorrelated Linear Discriminant Analysis (ULDA) to address these challenges, and offers a unified and stand-alone classifier. Through simulations and real-world datasets, the new framework demonstrates effective control of Type I error rates and improved classification accuracy, particularly in cases involving perfect group separations. The results highlight the potential of this approach as a robust alternative to the traditional stepwise LDA framework.




















[
    *
    Received: 20 June 2024 / Accepted: 05 July 2024
===================================================






§ INTRODUCTION


LDA seeks to find linear combinations of features that can best separate groups by maximizing the ratio of between-group variance to within-group variance. However, LDA is sensitive to noise variables and prone to overfitting. To address these issues, stepwise LDA is introduced, which iteratively adds or removes variables based on predefined inclusion and exclusion criteria. Various versions of stepwise LDA have been developed, ranging from stand-alone programs like DIRCRIM <cit.> and ALLOC-1 <cit.>, to options within statistical packages such as BMDP <cit.>, SPSS^® <cit.>, and SAS^® <cit.>. While specific implementations may differ in variable selection criteria, most follow a common framework discussed in <cit.>. Nonetheless, the heavy reliance on Wilks' Λ presents several challenges, some of which can be mitigated by substituting it with Pillai's trace.


Traditional LDA relies on the inverse of the within-class scatter matrix, leading to computational issues when the matrix is non-invertible. In contrast, ULDA <cit.> uses a different loss function based on the trace to solve this problem. Since both ULDA and Pillai's trace use trace-based criteria, it is logical to integrate them to develop a more robust stepwise LDA framework.


This paper is organized as follows: Section <ref> discusses the limitations of Wilks' lambda in the traditional stepwise LDA framework and highlights the computational and statistical challenges that arise. Section <ref> introduces our proposed forward selection framework based on Pillai's trace and ULDA, focusing on the algorithmic advancements and theoretical properties developed to resolve the challenges discussed. In Section <ref>, we present empirical analyses, including simulations and real data analyses, to demonstrate the effectiveness of the proposed method in controlling Type I error rates and improving classification accuracy when Wilks' Λ fails. We conclude in Section <ref>.







§ PROBLEMS WITH WILKS' Λ



First, we briefly introduce the most widely used stepwise LDA framework. Suppose we have a data matrix 𝐗∈ℝ^N × M with N observations and M features. Our response 𝐲∈ℝ^N is a factor vector containing J classes. Let 𝐱_ji∈ℝ^M represent the ith observation from class j, 𝐱̅_j ∈ℝ^M be the mean vector for class j derived from its n_j instances, and 𝐱̅∈ℝ^M denote the overall mean vector across all samples. 𝐇_B ∈ℝ^J × M and 𝐇_W ∈ℝ^N × M are defined as:


    𝐇_B    = [√(n_1)(𝐱̅_1-𝐱̅), √(n_2)(𝐱̅_2-𝐱̅), …, √(n_J)(𝐱̅_J-𝐱̅)]^T, 
    𝐇_W    = [(𝐱_11 - 𝐱̅_1), …, (𝐱_1n_1 - 𝐱̅_1), (𝐱_21 - 𝐱̅_2), …, (𝐱_2n_2 - 𝐱̅_1), …, (𝐱_Jn_J - 𝐱̅_J)]^T.



Then, the between-class scatter matrix 𝐒_B, within-class scatter matrix 𝐒_W, and total scatter matrix 𝐒_T can be defined as:


𝐒_B    = ∑_j=1^J n_j(𝐱̅_j-𝐱̅)(𝐱̅_j-𝐱̅)^'        = 𝐇_B^T𝐇_B 

𝐒_W    = ∑_j=1^J ∑_i=1^n_j(𝐱_ji-𝐱̅_j)(𝐱_ji-𝐱̅_j)^'        = 𝐇_W^T𝐇_W 

𝐒_T    = ∑_j=1^J ∑_i=1^n_j (𝐱_ji-𝐱̅)(𝐱_ji-𝐱̅)^'        = 𝐒_B + 𝐒_W.



Let 𝐒_T(1,2,…,p) and 𝐒_W(1,2,…,p) be the total and within-class scatter matrix with p variables {𝐱^(1), 𝐱^(2), …, 𝐱^(p)} added. Then the Wilks' Λ is defined as:


    Λ(1,2,…,p) = (𝐒_W(1,2,…,p))/(𝐒_T(1,2,…,p)).


After adding 𝐱^(p+1), we use partial Wilks' Λ to evaluate its marginal effect:


    Λ(p+1)=Λ(1,2, …, p, p+1)/Λ(1,2, …, p).


The null hypothesis H_0 states that the variables {𝐱^(1), 𝐱^(2), …, 𝐱^(p+1)} are from a multivariate normal distribution and are independent of the response 𝐲. Unless otherwise specified, this H_0 will be assumed as the null hypothesis throughout the remainder of this paper. Under H_0, the partial F-statistic follows an F-distribution:


    F=N-J-p/J-11-Λ(p+1)/Λ(p+1)∼ F_J-1, N-J-p.


In the (p+1)-th step, partial F-statistics are calculated for the remaining M-p variables, and the variable with the largest F-statistic is selected. It will be included in the model if it meets specific inclusion criteria, such as F ≥ 4, or if the corresponding p-value is below α.


Following the addition of a variable, the deletion phase begins. With p+1 variables now in the model, p+1 new pairs of scatter matrices (𝐒_W_i, 𝐒_T_i) are computed, each excluding one variable 𝐱^(i). The partial F-statistics are then calculated for each pair, and the variable associated with the smallest F-statistic is considered for removal if the exclusion criterion is satisfied (e.g., F < 3.996 in BMDP). This stepwise process continues until all variables have been added, or no further variables can be added or removed.


Next, we introduce three major drawbacks of using Wilks' lambda in the current stepwise LDA framework.





 §.§ Premature Stopping



When perfect linear dependency exists in the data matrix, we would expect 0/0 on the right-hand side of equation (<ref>), causing errors in some stepwise LDA programs, such as  in R. Wilks' Λ is not well-defined under perfect linear dependency, and to allow the stepwise selection to continue, a quick fix is to manually set it to 1, indicating no discrimination power.


We know from equation (<ref>) that the partial Λ is the ratio of two Wilks' Λ. Most programs will stop the stepwise LDA process when Λ = 0, as all subsequent partial Λ calculations become 0/0 and are thus ill-defined. However, stopping at Λ = 0 isn't always appropriate, as it indicates that one group of classes is perfectly separable from another, but it doesn't necessarily imply perfect separation of all classes in non-binary classifications, as shown in Figure <ref>. After selecting X_2, Wilks' Λ = 0 since the within-class variance is zero on X_2, causing the stepwise selection to stop. It successfully separates class A from classes B and C but cannot distinguish class B from class C. Additionally, when multiple variables result in Λ = 0, only one is selected, leading to the potential waste of useful information contained in the remaining variables.








 §.§ Partial F's Distribution Under Stepwise Selection



Here, we use an example to demonstrate that the distribution of the partial F-statistic does not follow an F-distribution under the stepwise selection framework, casting doubt on the validity of the associated hypothesis testing. Under H_0, the partial F-statistic follows an F-distribution <cit.>. However, the original proof assumes that variables are ordered randomly, rather than selected through a stepwise process. Intuitively, if we maximize the F-statistic at each step, the result will be biased, as noted in <cit.>.


The simulation setup is as follows: in each round, we simulate N = 150 observations from J = 3 classes, with each class having the same sampling probability of 1/3. We simulate X_1 from a standard normal distribution in the one-variable scenario and simulate X_1 and X_2 from independent standard normal distributions in the two-variable scenario. We simulate 10,000 rounds and record the partial F-statistic from each round at each step. We then compare the simulated F-statistic with the theoretical distribution, as summarized in Figure <ref>.


 The upper plot corresponds to the one-variable scenario. Since we have only one variable, X_1, no selection occurs, and the simulated distribution matches the theoretical distribution closely. The lower plot shows the two-variable scenario. Here, the stepwise selection first chooses X^(1), which has a larger partial F-statistic (and a smaller partial Λ-statistic) compared to the other variable, resulting in an upward bias in the first partial F-statistic. The second partial Λ is the ratio of the overall Wilks' Λ (with two variables) to the first (partial) Wilks' Λ, so its partial Λ-statistic is biased upwards and its partial F-statistic is biased downwards. Note that the theoretical distributions for the first and second partial F are different (F_2,147 and F_2,146), but the difference is negligible, so we assign the same color to both distributions in the plot.







 §.§ Inflated Type I Error Rate



In most programs, a fixed threshold of 4 is applied to the partial F-statistic. Another possible criterion is comparing the p-value of the partial F-statistic with the predefined α. Here, we simulate two scenarios to demonstrate that the type I error is inflated since both methods fail to account for the number of variables screened. For simplicity, forward selection is used instead of stepwise selection throughout this paper unless otherwise stated.


We use the iris dataset for our first simulation. It contains N = 150 flowers from J = 3 species (50 setosa, 50 versicolor, and 50 virginica), along with four features that characterize the flowers. We then add M mutually independent standard normal noise variables to it. Stepwise selections with both types of thresholds are performed, and we conclude that a type I error is made if at least one noise variable is selected. We let M = 1, 2, 4, 8, 16, 32, 64, 128, and for each M we repeat the simulation 2,000 times to obtain a confidence band. α is set to 0.05 throughout this paper unless otherwise specified.


In the second scenario, we simulate the null case where no variables are informative. We reuse the setup from the first scenario but remove all four original features from the iris dataset, leaving all remaining features independent of the species. The results from both scenarios are summarized in Figure <ref>. Both methods fail to control the type I error in either scenario. Due to the issue with multiple testing, the type I error rate increases as the number of noise variables grows. The fixed threshold of 4 performs slightly better than using the p-value, partly because a threshold of 4 corresponds to a p-value of approximately 0.02 in this setting, based on F_2,147(4) ≈ F_2,143(4) ≈ 0.02.







§ THE PROPOSED ALGORITHM



In this section, we first introduce ULDA and compare it to classical LDA. Next, we present our enhancements to ULDA and make it a stand-alone classifier. We then introduce the forward ULDA framework and derive the distributions of the test statistics used. Finally, we demonstrate that the type I error rate is well controlled within the new framework and summarize the algorithm.




 §.§ ULDA: Overview and Enhancements



Uncorrelated LDA (ULDA) is an extension of LDA that addresses scenarios where the within-class scatter matrix, 𝐒_W, is not invertible. Fisher's criterion aims to find transformation vectors 𝐰∈ℝ^M that maximizes the ratio:


    max _𝐰𝐰^T 𝐒_B 𝐰/𝐰^T 𝐒_W 𝐰



The optimal 𝐖 is derived by solving an eigenvalue decomposition on 𝐒_W^-1𝐒_B. The resulting eigen vectors 𝐖 = [𝐰_1, 𝐰_2, …] projects the original data 𝐗 into orthogonal linear discriminant scores 𝐗𝐰_i, ranked in descending order of their signal-to-noise ratios (eigenvalues). However, challenges arise when 𝐒_W is not invertible, such as when there are more variables than observations or when variables are linearly dependent. On the other hand, ULDA uses a different criterion:

    𝐖=max _𝐖^T 𝐒_T 𝐖=Itrace((𝐖^T 𝐒_T 𝐖)^+(𝐖^T 𝐒_B 𝐖)),


where 𝐀^+ denotes the Moore-Penrose inverse of matrix 𝐀. This ensures that the ULDA solution always exists, and <cit.> shows that ULDA is equivalent to classical LDA when 𝐒_T is nonsingular. Several insights can be drawn from equation (<ref>):


	
  * 𝐖^T 𝐒_T 𝐖 = I, meaning we deliberately discard the null space of 𝐒_T. This approach is reasonable because the null space of 𝐒_T is the intersection of the null spaces of 𝐒_W and 𝐒_B, and equation (<ref>) is not affected by vectors from the null space of 𝐒_B.
	
  * If we change the constraint from 𝐖^T 𝐒_T 𝐖=I to 𝐖^T𝐖=I, it becomes Orthogonal LDA (OLDA) <cit.>. However, the constraint 𝐖^T 𝐒_T 𝐖=I is more beneficial when constructing an LDA classifier. To solve equation (<ref>), <cit.> presents an algorithm based on Generalized Singular Value Decomposition (GSVD), where 𝐒_B, 𝐒_W, and 𝐒_T are diagonalized simultaneously. Suppose the rank of 𝐒_T is M, then 


    𝐖^T 𝐒_B 𝐖   = diag(α_1^2, α_2^2, …, α_M^2) 
    𝐖^T 𝐒_W 𝐖   = diag(β_1^2, β_2^2, …, β_M^2) 
    𝐖^T 𝐒_T 𝐖   = diag(α_1^2 + β_1^2, α_2^2 + β_2^2, …, α_M^2 + β_M^2) 
       = diag(1, 1, …, 1)

	
	Since the likelihood-based LDA classifier depends on the inverse of 𝐒_W, computational resources are saved if it has already been diagonalized.
	
  * Some commonly used test statistics from MANOVA are related to equation (<ref>). Pillai's trace is defined as V = trace(𝐒_T^-1𝐒_B), which equals ∑_i=1^M α_i^2 after transformation. In other words, ULDA can be viewed as maximizing the generalized Pillai's trace under certain constraints. On the other hand, Wilks' Λ = ∏_i=1^M β_i^2. We discussed the issue with Λ = 0 in Section <ref>. Λ = 0 means that β_i = 0 and α_i = 1 for some i. In the stepwise selection framework, since Wilks' Λ is a product, once it becomes zero, it remains zero. On the other hand, Pillai's trace is a summation, and adding another 1 has no side effect. Pillai's trace is also superior to Wilks' Λ in other aspects <cit.>.


Next, we introduce our speed enhancement for the ULDA algorithm when N > M. <cit.> presents a ULDA algorithm that diagonalizes 𝐒_T and 𝐒_B separately. Based on our experience, it is slower by a constant factor compared to the GSVD-based version <cit.>, which is described in Algorithm <ref> (rewritten to suit our needs). However, when the sample size N is large, the SVD decomposition on 𝐊∈ℝ^(J+N) × M in the line <ref> of Algorithm <ref> creates a runtime bottleneck. This can be resolved by reducing the dimension of 𝐊 before performing SVD (or complete orthogonal decomposition). Since 𝐇_W contributes most of the dimensionality, and the SVD depends on 𝐇_W^T𝐇_W, one possible solution is to replace 𝐇_W with 𝐆_W ∈ℝ^M × M, where 𝐇_W^T𝐇_W = 𝐒_W = 𝐆_W^T𝐆_W. We suggest performing a reduced QR decomposition 𝐇_W = 𝐐_W𝐑_W and replacing 𝐇_W with 𝐑_W, so that we have

    𝐇_W^T𝐇_W = 𝐑_W^T𝐐_W^T𝐐_W𝐑_W = 𝐑_W^T𝐑_W.


<cit.> follows a similar approach, using a Cholesky decomposition 𝐒_W = 𝐂_W^T𝐂_W and replacing 𝐇_W with 𝐂_W. We now use a simulation to evaluate the performance of these two variants and the original ULDA.




The simulation setup is as follows: in each round, we simulate N = 10000 observations from J = 10 classes, with each class having the same sampling probability of 1/10. The features are M mutually independent standard normal noise variables. We then use the three algorithms to calculate the transformation matrix 𝐖. We let M = 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, and for each M we repeat the simulation 30 times to obtain a confidence band. The results are summarized in Figure <ref>. Their differences in runtime become larger as the number of features increases. For a data matrix of dimension 10000 × 1024, ULDA with QR decomposition is 38.6% faster than the original GSVD-based ULDA implementation (6.2 seconds vs. 10.1 seconds). Consequently, we added this additional QR decomposition step into the ULDA pipeline.




Another enhancement is our integration of the likelihood structure. The original ULDA method primarily serves as a dimension reduction tool and leaves the classification task to K-nearest-neighbors or other classifiers. However, LDA assumes Gaussian density for the features, and posterior probabilities naturally serve as a powerful classifier, as implemented in the R package  and the Python package . One difficulty lies in the invertibility of 𝐒_W. Occasionally, we find a highly discriminative direction such that α^2 = 1 (or Wilks' Λ = 0), meaning the ratio of the between-class variance to the total variance is 1, and one group of classes is perfectly separable from another. We address this by manually setting the within-class variance to 10^-5, which has proven effective in our experience. To understand why it works, note that the discriminant function in LDA is:


    δ_j(𝐱)=𝐱^T Σ_W^-1μ_j-1/2μ_j^T Σ_W^-1μ_j+logπ_j,


where π_j is the prior for class j, and we have μ_j = 𝐱̅_j and Σ_W = 𝐒_W / (N - J). 𝐱 will be classified to class k if k = max_j δ_j(𝐱). When perfect separation occurs, we aim to capture this through the discriminant function. By setting the within-class variance to 10^-5, the Mahalanobis norm associated with that direction is magnified 10^5 times, which is large enough to dominate the effects from all other directions.


For missing values, we impute them with the median for numerical variables and assign a new level for categorical variables. Additionally, we include missing value indicators for numerical variables. For categorical variables, we use one-hot encoding and transform them into dummy variables. We also accommodate unequal misclassification costs. Let 𝐂 represent the misclassification costs, where C_ij = C(i | j) is the cost of classifying an observation into class i given that it belongs to class j. Suppose the predicted posterior probability for an observation is 𝐩 = (p̂_1, p̂_2, …, p̂_j); then the cost of predicting it to class i is


    C_i = ∑_j = 1^J C(i | j) * p̂_j,


and 𝐱 will be classified to class k if k = max_i C_i(𝐱).





 §.§ Forward ULDA: Distribution and Threshold



We now derive the distribution of the test statistics used in the forward ULDA framework. Without loss of generality, we assume 𝐒_T is invertible throughout this section, as redundant columns that cause 𝐒_T to be non-invertible have no discriminative power and can always be removed.



Pillai's trace is non-decreasing when new variables are added to the model.



Suppose 𝐗∈ℝ^N × K have been included in the model, and the new variable to be added is 𝐳∈ℝ^N. The new between-class and total scatter matrices with (K+1) variables can be written as:


    𝐒_B =
    [   𝐁_x   𝐛_x; 𝐛_x^'   b_z ]  𝐒_T =
    [   𝐓_x   𝐭_x; 𝐭_x^'   t_z ],


where 𝐁_x and 𝐓_x are the previous between-class and total scatter matrices for 𝐗. If K = 0, the difference in Pillai's trace will be b_z/t_z. This value is non-negative since the between-class scatter matrix is positive semi-definite (b_z ≥ 0) and the total scatter matrix is positive definite (t_z > 0). For K ≥ 1, we aim to show


    trace(𝐒_T^-1𝐒_B) - trace(𝐓_x^-1𝐁_x) ≥ 0.


Since 𝐒_T and 𝐓_x are invertible, we have











    !𝐒_T^-1 =
    [ 𝐓_x^-1 + 𝐓_x^-1𝐭_x𝐭_x^'𝐓_x^-1(t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1                    -𝐓_x^-1𝐭_x(t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1;                  -(t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1𝐭_x^'𝐓_x^-1                              (t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1 ].


According to the block matrix multiplication and the properties of the trace, we have


    trace(𝐒_T^-1𝐒_B)    = trace(𝐓_x^-1𝐁_x) 
          + (𝐭_x^'𝐓_x^-1, -1)𝐒_B(𝐭_x^'𝐓_x^-1, -1)^'
         × (t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1.


(t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1 is the Schur complement of the block 𝐓_x of the matrix 𝐒_T. Since 𝐓_x and 𝐒_T are both positive definite, we have (t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1 > 0. (𝐭_x^'𝐓_x^-1, -1)𝐒_B(𝐭_x^'𝐓_x^-1, -1)^' is a quadratic form, and since the middle matrix is positive semi-definite, it is non-negative.


At its core, ULDA seeks to maximize Pillai's trace V = trace(𝐒_T^-1𝐒_B). According to Theorem <ref>, with each variable added, the current Pillai's trace increases (or remains the same). Let V^(M) denote the Pillai's trace with all M variables included. The goal of forward selection is to approximate V^(M) using V^(K), where K ≪ M.


Suppose the variable set {𝐱_1, 𝐱_2, …, 𝐱_K-1} has been selected after the first (K-1) steps, and the Pillai's trace of that variable set is V^(K-1)_max. Here, the subscript indicates that this Pillai's trace is not of (K-1) randomly selected variables but is instead maximized at each step through forward selection. At step K, we calculate V^(K)_(1), V^(K)_2, …, V^(K)_(M-K+1), where V^(K)_(i) denotes the Pillai's trace of the variable set {𝐱_1, 𝐱_2, …, 𝐱_K-1, 𝐱_(i)}. Let k = max_i V^(K)_(i). We then select 𝐱_(k) as the best candidate at step K, and V^(K)_max = V^(K)_(k). To establish an inclusion criterion, we must measure the marginal effect of the added variable 𝐱_(k), which corresponds to V^(K)_max - V^(K-1)_max.



At step K, let t_K be the (1 - α)^1/(M-K+1) quantile of B(J-1/2, N-J/2). Then, P(V^(K)_max - V^(K-1)_max≥ t_K) ≤α as N →∞ under H_0, where the newly added variable 𝐳 is normally distributed and independent of both 𝐗 and 𝐲.



When K = 1, V^(1) under H_0 follows a beta distribution B(J-1/2, N-J/2) <cit.>. Since V^(1)_max - V^(0)_max = V^(1)_max is the maximum of V^(1)_1, V^(1)_2, …, V^(1)_M, which are M i.i.d. random variables from the beta distribution, its CDF can be written as I_x^M(J-1/2, N-J/2) where I_x(J-1/2, N-J/2) is the CDF of B(J-1/2, N-J/2). To control the type I error below α, the threshold t must satisfy I_t^M(J-1/2, N-J/2) ≤ 1 - α, which is equivalent to I_t(J-1/2, N-J/2) ≤ (1 - α)^1/M. Then t is the (1 - α)^1/M quantile of B(J-1/2, N-J/2). 


When K > 1, from equation (<ref>) we know that


    V^(K)_(i) - V^(K-1)_max   = (𝐭_x^'𝐓_x^-1, -1)𝐒_B(𝐭_x^'𝐓_x^-1, -1)^'
         × (t_z - 𝐭_x^'𝐓_x^-1𝐭_x)^-1.


This equation still holds if we replace 𝐒_B and 𝐒_T with 𝐒_B/(N-J) and 𝐒_T/(N-J), which are the least squares estimators of the between-class and total covariance matrices. Since 𝐗 and 𝐳 are independent, their covariance 𝐭_x →0 as N →∞. Note that 𝐒_B/(N-J) and 𝐒_T/(N-J) are finite as N →∞. Substituting 𝐭_x = 0 into equation (<ref>), we get


    V^(K)_(i) - V^(K-1)_max   = (0, -1)𝐒_B(0, -1)^'(t_z - 0)^-1
       = b_z/t_z


b_z/t_z is Pillai's trace for 𝐳. Therefore, the distribution of V^(K)_(i) - V^(K-1)_max can be approximated by V^(1), and the rest follows the scenario where K = 1.


Based on our experience, this asymptotic approximation sometimes leads to a very conservative threshold, with the type I error falling well below the predefined α. Therefore, we introduce an empirical approximation to mitigate this problem and achieve higher power. Suppose we have already added K-1 variables and the current Pillai's trace is V^(K-1)_max. Since the maximum Pillai's trace for J classes is J-1, the maximum Pillai's trace that can be added is bounded by J - 1 - V^(K-1)_max, which can be viewed as the maximum Pillai's trace for a classification problem with J - V^(K-1)_max classes. Thus, at the k-th step, the threshold becomes the quantile from B(J^'-1/2, N-J^'/2) instead of B(J-1/2, N-J/2), where J^' = J - V^(K-1)_max.





 §.§ Type I Error: Analysis and Control



Here, we analyze the type I error under the forward ULDA framework and demonstrate that the family-wise error rate is controlled at the nominal level α. Suppose we have M variables in total, some of which are noise variables (𝐱∈ S_n) and some are informative (𝐱∈ S_i). At each step, there are three possible outcomes: a noise variable is selected, the selection stops, or an informative variable is selected. The entire process is illustrated in Figure <ref>.




Suppose 𝐱 is the variable with the largest Pillai's trace and is selected at the K-th step. Conditional on whether 𝐱∈ S_n or 𝐱∈ S_i, there are four possible outcomes:


    p_K1   = P(𝐱 is added|𝐱∈ S_n) 
    
    	p_K2   = P(𝐱 is not added|𝐱∈ S_n) 
    
    	p_K3   = P(𝐱 is added|𝐱∈ S_i) 
    
    	p_K4   = P(𝐱 is not added|𝐱∈ S_i)


In situations with p_K2 and p_K4, 𝐱 is not added, and forward selection stops. Therefore, no type I error is made or will be made. This corresponds to the green regions in Figure <ref>. For the situation with p_K3, since an informative variable is added, no type I error is made at the current step, corresponding to the yellow regions in Figure <ref>. The purple region in Figure <ref> reflects scenarios where a type I error is made, with p_K1 being the only situation that results in such an error. Theorem <ref> shows that under H_0, p_K1≤α, meaning that at each step, the probability of branching into the purple region is controlled at α. Now, we aim to show that, overall, the probability of ending up in any purple region is controlled at α.


The probability of reaching node 2 is p_11≤α. The probability of reaching node 5 is p_13× p_21≤ p_21≤α. The reason we can use the product of p_13 and p_21 to calculate this probability is that under H_0, the variable selected in the first step is assumed to be independent of the variable selected in the second step. For nodes like node 2 and node 5, where the first noise variable is added in the current step, the probability of reaching them can be written as


    p_K1∏_k = 1^K-1p_k3≤ p_K1≤α.


Meanwhile, the probability of reaching their child nodes is also controlled at α, because reaching these nodes requires first reaching their parent node. All purple nodes fall into one of these two scenarios, so the family-wise type I error rate is controlled at α. This means that if the forward ULDA selects K variables {𝐱_(1), 𝐱_(1), …, 𝐱_(K)}, then the probability that at least one 𝐱_(i) is a noise variable is controlled at α.


The forward selection framework is summarized in Algorithm <ref>.








§ EMPIRICAL ANALYSIS



In this section, we use simulation and real data to showcase the performance of three forward LDA variants:


	
  * : the proposed variant using Pillai's trace (see Algorithm <ref>). 
	
  * : the original variant based on Wilks' Λ. The inclusion criterion is based on p-value, with a variable included if the p-value is below the predefined α. 
	
  * : This variant applies an additional Bonferroni correction to the p-value compared to the second variant. If there are (M-K+1) variables to choose from at the K-th step, the p-value is multiplied by (M-K+1) to adjust for the multiple testing.


Note that these variants are for selection alone. They help check the type I error and power (whether the desired variables are included). To further compare the testing accuracy, we apply ULDA (see Section <ref>) to the selected variables.





 §.§ Type I Error Evaluation on Iris: Pure Noise and Mixed Cases



We use the same simulation settings from Section <ref> to compare the three forward LDA variants, and the results are summarized in Figure <ref>.  and  successfully control the type I error in both scenarios. In contrast,  suffers from an inflated type I error rate due to multiple testing. These results validate Theorem <ref>, demonstrating that the type I error rate is well-controlled under H_0, where the noise variables are normally distributed and independent of both the informative variables and the response.








 §.§ Handling Λ=0: Analysis on Simulated and Real Data



In this section, we illustrate the primary advantage of our proposed method over the original framework: its ability to handle scenarios where Wilks' Λ = 0.


First, we use a simulated dataset, which contains 2,000 observations. The response variable is randomly selected from 10 classes, each with an equal probability of 1/10. We then create a dummy matrix (one-hot encoding) of the response, resulting in 10 columns, each consisting of 1s or 0s. These 10 columns are used as our features. Ideally, these 10 features can perfectly predict the response, and a robust forward selection method should select them all. However,  stops after selecting only one feature, I_Class One, the indicator of the first class. Here, I_Class One = 1 for observations from class one and I_Class One = 0 for the other classes. Therefore, the within-class variance is 0, leading to Wilks' Λ = 0. Meanwhile, with , the feature I_Class One contributes a value of 1 to the overall Pillai's trace, which does not trigger a stop.  continues adding features, with each feature I_Class i contributing a value of 1 to the overall Pillai's trace. It selects 9 features and then stops, as the maximum Pillai's trace of J - 1 = 9 is reached. With these 9 features, we have enough information to perfectly predict the response.


Next, we use real data from the National Highway Traffic Safety Administration (NHTSA) and its Vehicle Crash Test Database. This database contains data from various crash tests, including those conducted for research, the New Car Assessment Program (NCAP), and compliance purposes. The dataset captures a wide range of vehicle attributes, crash conditions, and safety outcomes. Our focus here is to predict the type of engine used in the tested vehicles. This dataset is challenging to analyze for several reasons:


	
  * Cyclic variables: Variables like impact angle are cyclic, where 359^∘ and 1^∘ should be very similar in real life but are very different in their numerical representation. We address this by transforming all angles to their cos and sin values.
	
  * Missing values: For missing values, we impute them with the median for numerical variables and assign a new level for categorical variables. Additionally, we include missing value indicators for numerical variables.
	
  * Multicollinearity: Variables BX1 to BX21 are measurements of different parts of the car, and some of them are highly correlated, which affects the performance of classical LDA, making forward ULDA more promising.


After data preprocessing, there are 3,273 crash tests and 173 (982 if all categorical variables are transformed into dummy variables) predictor variables. Our response variable, engine type, has 18 different types, with an imbalanced distribution: the most frequent class, 4CTF, has 1,250 occurrences, while the least frequent class, NAPP, has only one occurrence. We then apply both forward LDA variants,  and , to this dataset.


 stops after selecting variables I_Model=RX and I_Model=RX-8. The reason is that only these two Mazda models (RX and RX-8) use the ROTR engine type. Therefore, the within-class variance becomes zero, allowing it to perfectly identify the ROTR engine. However, this premature stopping prevents it from identifying all the other engine types. The order of the variables selected also reveals its lack of statistical power. The first selected variable should be the most discriminative. However, the top two variables selected by  can only correctly identify 3 instances out of 3,273, clearly not the most discriminative variables. This bias towards the perfect separation of certain classes, regardless of how few the instances are, dominates the forward selection process via smaller Wilks Λ, indicating that Wilks Λ is not an ideal test statistic and fails to capture the most discriminative information across all classes.


In contrast,  ends up selecting 175 variables. Using ULDA as the classifier and a 10-fold cross-validation, the prediction accuracies of  and  are 0.38 and 0.65, respectively. Results from simulation and real data demonstrate that  outperforms  when Wilks' Λ = 0.





§ CONCLUSION



In this paper, we present a new forward discriminant analysis framework based on Pillai's trace, which demonstrates superiority in situations involving perfect group separations (Wilks' Λ = 0) compared to traditional methods. Our approach effectively controls the type I error rate and integrates seamlessly with ULDA, providing a unified classifier.


Despite these advancements, there are some limitations to consider. The criterion used in forward selection, whether based on Wilks' Λ or Pillai's trace, serves only as a measure of goodness of fit within the MANOVA framework, as mentioned in <cit.>. It does not directly reflect the model's training or testing accuracy, although these two typically align well. For those primarily concerned with accuracy, we recommend using the forward selection framework to rank variables, followed by methods such as cross-validation to choose the best subset of variables.


Additionally, our derivation of the distribution of the test statistic in Section <ref> is asymptotic and includes an empirical approximation. There is potential to develop a more precise distribution for finite sample size under alternative assumptions, which could be a direction for future research.


In this paper, we focus primarily on forward selection rather than stepwise selection. The main reason is the difficulty in theoretically justifying a valid exclusion criterion. Moreover, based on our experience, the addition of a variable deletion step typically results in minimal improvement in classification accuracy while significantly increasing runtime.


The related R package  is available on CRAN.











Acknowledgements

The author gratefully acknowledges Prof. Wei-Yin Loh from UW-Madison for his unwavering guidance throughout the author's PhD journey.







§ DECLARATIONS


Conflict of interest: The author declares no conflict of interest.



































































plainnat





