The Linear Discriminant Analysis (LDA) stands out as a popular technique for classification tasks. However, its applicability falters when the sample is high-dimensional; that is, the number of input variables p is larger than the sample size n, which is a situation that has become increasingly prevalent in practical scenarios. In these cases, the classical LDA faces a critical limitation in estimating the population covariance matrix due to the singularity of the pooled sample covariance matrix. For example, Bickel and Levina (2004) pointed out that the performance of the classical LDA can be as poor as a random guess for a high-dimensional sample.

Many authors have studied modifications of the classical LDA concerning high-dimensional sample, with early proposals based on the independence rule and variable selection. Bickel and Levina (2004) proposed to use the independence rule or the diagonal LDA paradigm for estimating the population covariance matrix by ignoring the correlation structure of the input variables. Tibshirani et al. (2002) and Fan and Fan (2008) introduced the nearest shrunken centroid estimation and the features annealed independent rule, in which the variable selection is done by soft and hard thresholding rules, respectively. However, these methods do not achieve Bayes classification error for the classification task unless the true common covariance matrix is diagonal.

Recently, the sparse LDA approaches have been proposed as an alternative to the independence rule. Cai and Liu (2011) proposed the linear programming discriminant rule that finds a sparse estimate of the population covariance matrix by using the Dantzig selector (James et al., 2009), and Fan et al. (2012) proposed the regularized optimal affine discriminant method, which is based on the ℓ₁-Fisher’s discriminant analysis. Similar ℓ₁-penalized linear discriminant analysis was studied in Trendafilov and Jolliffe (2007) and Witten and Tibshirani (2011). Clemmensen et al. (2011) proposed the sparse optimal scoring method which solves the optimal scoring formulation (Hastie et al., 1994) with the ℓ₁-penalty. The Direct Sparse Discriminant Analysis (DSDA) studied in Mai et al. (2012) is another popular sparse LDA, which is computationally efficient and easier to understand, because it reformulates the high-dimensional LDA into a penalized linear regression framework. In general, the sparse LDA approaches have theoretical advantages over other methods of independence rules since they allow non-diagonal population covariance matrices, enabling the resulting classifier to achieve Bayes classification error. However, most previous works have been limited to specific convex penalties such as the Least Absolute Selection and Shrinkage Operator (LASSO) and its variants.

In this study, we focus on the DSDA applied with a class of concave penalties, including the smoothly clipped absolute deviation (Fan and Li, 2001), minimax concave (Zhang, 2010), and truncate d-ℓ₁ (Shen et al., 2013) penalties as examples. Within this framework, implementing the concave penalized direction vector of discrimination appears to be straightforward, but developing the theoretical properties of the estimated direction vector still remains challenging. For this purpose, we prove that the concave penalized direction vector satisfies an oracle property uniformly in the class, which is the main contribution of the paper. Here, the oracle property implies that an ideal LDA direction vector of discrimination can be exactly recovered through concave penalized least square estimation. In addition, we provide various numerical studies to confirm whether the theoretical results hold with finite samples.

The rest of the paper is organized as follows. Section 2 introduces the sparse LDA. Section 3 introduces the concave penalized LDA and presents the related theoretical results. Section 4 provides numerical studies to confirm the theoretical results, and concluding remarks are given in Section 5. Technical details and proofs are provided in Appendix.

2.1. Fisher’s linear discriminant analysis

Fisher’s Linear Discriminant Analysis (LDA) (Fisher, 1936) is an efficient technique for discriminating a binary class label C ∈ {0, 1} given an input vector X ∈ ℝ^p. The LDA assumes that X is a random vector with multivariate normal distribution given a binary random variable C such that

independently, where μ_c and ∑_c are the mean vector and covariance matrix of the normal distribution, respectively. Let φ(·; μ_c, ∑_c), c ∈ {0, 1} be the density function of N(μ_c, ∑_c) and π_c = P(C = c). The ratio between two conditional density values of C|X = x provides a natural rule for discriminating the class label as C = 1 given X = x as follows.

In addition, when the covariance matrices in (2.1) are assumed to be common, ∑₀ = ∑₁ = ∑, the discrimination rule in (2.2) can be simplified into a linear form,

where β^BS = ∑⁻¹(μ₁ – μ₀) is Bayes direction vector.

Let (x_i, c_i), i ≤ n be n random samples of (X, C) then the LDA estimate ψ̂ of ψ in (2.3) can be obtained by using simple moment type estimates:

where β̂^LDA = ∑̂⁻¹(μ̂₁ – μ̂₀) is the LDA direction vector,

and n_c is the number of samples with c_i = c.

2.2. Sparse linear discriminant analysis

In many real fields, application of the LDA raises two challenging problems because of sparsity of the model and high-dimensionality of the samples. The sparsity of the model implies that Bayes direction vector satisfies βjBS=0 for some j ≤ p but, in general, the LDA estimate in (2.4) produces β^jLDA≠0 for all j ≤ p. The high-dimensionality of the samples implies p > n, where the pooled sample covariance matrix ∑̂ in (2.5) becomes singular. In both cases, the discriminant analysis may produce unexpected bad results in prediction as well as interpretation. There are many literatures that figure out these problems and the penalized approach via the least squares estimation can be a nice solution (Mai et al., 2012).

Let y_i = (−1)^1+c_in/n_ci, i ≤ n, c_i ∈ {0, 1} then the LDA direction vector in (2.4) can be interpreted as a Least Squares Estimator (LSE) as follows:

for some positive constant k unless μ̂₀ ≠ μ̂₁ (Hastie et al., 2009), where

Hence the LDA rule in (2.4) can be cast into the LSE rule:

Given the sparsity and high-dimensionality, it is natural to employ the penalized LSE:

where J_λ is a penalty equipped with a tuning parameter λ. In this case, the LSE rule in (2.8) can be replaced with the penalized LSE rule again:

where k_λ = β̂^λT∑̂ β̂^λ/(μ̂₁ – μ̂₀)^Tβ̂^λ is given by Mai et al. (2012) that is optimal to the discrimination whenever ( μ̂₁ – μ̂₀)^Tβ̂^λ > 0.

3.1. Concave penalized estimation

Let J_λ be a penalty function with tuning parameter λ > 0 and Jλ′ be the first derivative of J_λ. We assume that following conditions hold.

• (J1) J_λ(t) is concave and non-decreasing over t ∈ [0, ∞) and J_λ(0) = 0.

• (J2) Jλ′(t) is non-increasing and continuous over t ∈ (0, ∞) and limt→0+ Jλ′(t)=λ.

• (J3) Jλ′(t)≥λ-t/a over t ∈ (0, aλ) and Jλ′(t)=0 over t ∈ (aλ, ∞) for some a > 0.

The class of penalties that satisfy (J1), (J2), and (J3) has been studied as a representative concave penalty class (Kim and Kwon, 2012; Zhang and Zhang, 2012), including the Smoothly Clipped Absolute Deviation (SCAD) (Fan and Li, 2001),

minimax concave (MCP) (Zhang, 2010),

and truncated-ℓ₁ penalties (TLP) (Shen et al., 2013)

as examples, where x₊ = max{0, x}.

Let Z = (Z₁, . . . , Z_p) be the centered design matrix in (2.7) then the LSE direction vector in (2.7) can be simplified into

Now, considering the true signal set S={j∣βjBS≠0}, one of the best estimators for the sparse Bayes direction vector is the oracle LSE,

which is unavailable in practice without the knowledge of . Hence, the main objective of the concave penalized LDA is to recover β̂^OR through the penalized LSE:

where

3.2. Optimality conditions

In this subsection, we introduce three lemmas that provide non-asymptotic sufficient conditions for a given estimator to be a local or unique local minimizer of L_λ. Let Ξ_λ be the set of all local minimizers of L_λ and θ̂ = μ̂₁ – μ̂₀.

3.3. Oracle property

In this subsection, we present the main results of the paper: The concave penalized LSE in (3.2) is unique, and hence, the global minimizer of L_λ, and the same as the oracle LSE in (3.1) asymptotically, in the sense that Ξλκ={β^λ} and β̂^λ = β̂^OR with probability tending to 1.

3.3.1. Notations

We introduce some notations. For any matrix A, let ||A||_∞ = max_i ∑_j |A_ij|, where A_ij is the (i, j) entry of A, and let be a sub-matrix of A constructed by the entry A_ij, . For any vector a, let ||a||_∞ = max_i |a_i| and ||a||₁ = ∑_i |a_i|, where a_i is the ith element of a, and let be a sub-vector of a constructed by the elements a_i, . Let be the cardinality of a set . For any positive sequences x_n and y_n, we write x_n ≫ y_n if x_n/y_n → ∞ and x_n ≍ y_n if x_n/y_n → a for some constant a as n→ ∞.

Recall the definition of the oracle LSE β̂^OR obtained under the knowledge of the true signal set S={j∣βjBS≠0}. It is easy to see that

where Ω̂ = Z^TZ/n, θ̂ = μ̂₁ – μ̂₀, and N={j∣βjBS=0} is the true noisy set. Hence, we can construct a population counterpart β^OR of β̂^OR by defining

where Ω = Cov(X) and θ = μ₁ – μ₀. Note that β^OR is exclusively considered for the development of theoretical properties since β^OR = kβ^BS for some constant k ≠ 0. For further details, refer to Proposition 3 in Mai et al. (2012).

3.3.2. Oracle property

The main objective of the concave penalized LSE is to recover β̂^OR by using β̂^λ in (3.2), which is often referred to the strong oracle property (Fan and Li, 2001; Kim et al., 2016). From the lemmas in previous subsection, it is sufficient to prove that β̂^OR satisfies the outlined uniqueness and Sparse Riesz conditions with probability tending to 1 of which the proofs are provided in Appendix. We assume the following regularity conditions:

(C1) There exist positive constants, b_i, i ≤ 4, such that

for any n.

In this section, we report some results of simulation studies and real data analysis.

4.1. Simulation studies

We generated n random samples of C from Ber(π) with π = 1/2, and then generated n random samples of X|C = c, c ∈ {0, 1} from N(μ_c, ∑). We set μ₀ = 0 and μ₁ = ∑β^BS, where ∑_jk = r^|j–k| and βjBS=(2/j)(-1)jI(j≤q), j, k ≤ p. We consider three concave penalties, TLP with a = 0.001, MCP with a = 1.001, and SCAD with a = 2.001, as examples in the penalty class, and set n ∈ {500, 1000, 2000}, p ∈ {300, 3000}, and q ∈ {10, 20} with r ∈ {0.25, 0.5} for each simulation. We repeated the simulation 200 times by using R package ncpen (Kim et al., 2020).

We first investigated whether the oracle property can hold with finite samples and Table 1 shows the estimated probabilities of achieving the oracle property. For each simulation, we first found the interval [λ_max, λ_min] that satisfies ||β̂^λ_max||₀ = 0 and ||β̂^λ_min||₀ = 50. Then we checked whether there exists a λ that satisfies β̂^λ = β̂^OR by investigating 200 values of λ decreasing with log-scale in the interval. From the table, we observed the followings:

• The ratios of the TLP and MCP approaches nearly 1 for some cases, while the SCAD has lower ratios, with the largest ratio being 0.79. In cases where both p and q are large, it is rare for the oracle property to hold, but we believe that this occurs due to limitations in the simulation settings. In general, we can conclude that the ratio increases as n increases and p, q, and r decreases, regardless of the simulation settings and penalties, which confirms the result of Theorem 1.

• In addition, we checked the sign consistency, sign(β̂^λ) = sign(β̂^OR) for some λ, which is a slightly less stringent condition than the oracle property. The results exhibit similar pattern to that of the oracle property but the ratios are slightly larger. We found that the LASSO also achieves the sign consistency as Mai et al. (2012) proved. In general, the sign consistency of the LASSO does not hold even for the ordinary linear regression since it requires the strong irrepresentale condition on the design matrix. See Section 3 and example 3 in Zhao and Yu (2006) for some details.

Second, we compared the finite sample performance of the concave penalized estimators, using the LASSO as a benchmark method. The primary objective of the simulation is to check whether the oracle property can be realized through typical tuning parameter selection criteria. Note that there are two natural tuning parameter selection criteria from the characteristic of the framework: Minimizing the regression error in (2.7) and the classification error in (2.9). We used n training samples for estimating the penalized direction vector and n/2 independent validation samples for selecting an optimal tuning parameter λ_opt by minimizing the two criteria. We calculated three measures of selection accuracy: The numbers of true and false positive selection (TPS and FPS) and the ratio of the correct model selection (CMS): TPS=ΣjI(β^jλopt≠0,βjBS≠0), FPS=ΣjI(β^jλopt≠0,βjBS=0), and CMS is the ratio of the cases when TPS = q and FPS = p–q. In addition, we calculated the miss-classification error rate (ERR) by using 2n independent test samples.

Tables 2 and 3 summarizes the results. We first discuss the cases when λ_opt was chosen based on the regression errors of the validation samples:

• TPS increases to q as n increases, regardless of the penalties and simulation settings. This suggests that the signal variables tend to be selected as the sample size increases. As p or r increases while keeping n fixed, TPS decreases. This implies that higher model complexity and stronger correlation among the input variables have a negative effect on selection of signal variables. Among the concave penalties, the SCAD produces the largest TPS but the difference does not seem significantly large.

• For the concave penalties, FPS decreases as n increases for almost all cases. This implies that the noisy variables tend to be excluded as the sample size increases, which is not true for the LASSO. FPS of the SCAD is smaller than that of the LASSO but it is not at a negligible level even when n = 2000. However, FPS for the TLP and MCP are near 0 when n = 2000, implying that they excluded nearly all noisy variables.

• For the concave penalties, CMS increases as n increases regardless of the penalties and simulation settings, achieving the largest score 0.80, 0.77, and 0.34 for the TLP, MCP, and SCAD, respectively. For the LASSO, CMS is consistently 0, indicating that the regression errors does not work for the LASSO to select correct models.

• One of the reasons for the poor FPS and CMS of the SCAD seems to be because its shape or concavity on the interval [0, aλ]. Since the SCAD penalty is the same as the LASSO on [0, λ] and a > 2 is limited below, the maximum concavity on [0, aλ] is smaller than that of the MCP with a > 1, which results in selecting more variables with higher FPS and smaller CMS than those of the MCP. We refer to Zhang (2010) for some detailed discussion on the maximum concavity of the concave penalties.

• ERR decreases to that of the Bayes, the best performance, as n increases but increases as p or r increases with fixed n, regardless of the penalties and simulation settings.

For the cases of classification errors, we observed the followings:

• TPS shows similar patterns as n increases, resembling the case of regression errors. This occurs regardless of the penalties and simulation settings.

• The LASSO and SCAD select significantly less noisy variables, producing lower FPS and slightly higher CMS. However, the TLP and MCP show opposite results.

To sum up, (a) we can conclude that tuning parameter selection based on the regression errors exhibits better selection performance for the TLP and MCP. However, for the LASSO and SCAD, the classification errors seems to be more informative, as indicated by smaller FPS and larger CMS. (b) The two tuning parameter selection criteria can function to some extent in the simulation settings designed in this paper. (c) However, the probabilities of correctly identifying the true model in Table 1 consistently reach nearly 1, regardless of the penalties and simulation settings. This suggests that we need to develop other alternatives such as the information criterion based on the underling probability structure.

4.2. Real data analysis

We apply the penalized LDA methods to two benchmark datasets: The prostate and lung cancer datasets, which is available from the R packages datamicroarray. The prostate cancer dataset consists of the expression levels of approximately 12,600 genes obtained from 52 prostate tumour samples and 50 non-tumour prostate samples (Singh et al., 2002). Likewise, the lung cancer dataset comprises 12,533 genes from 150 patients with adenocarcinoma and 31 patients with malignant pleural mesothelioma (Gordon et al., 2002). The main task is to predict whether an observation is the specific tumor tissue or not, and to identify the genes associated with each type of cancer. We first choose top p ∈ {500, 1000, 2000} genes with the largest t-statistics across two classes. We then applied the penalized LDA methods with two classes as a response variable and the chosen top p genes as predictive variables. For comparison, we split the data into the training and test sets by randomly choosing 3/4 samples and 1/4 samples, respectively. For each training set, the tuning parameter λ is selected by the 10-fold cross-validation methods based on regression error and classification error as in the simulations. We compute the mis-classification error rate (ERR), the number of mis-classified samples (#MIS) on the test set, and the model size (SIZE) representing the number of selected variables on the training set.

Table 4 presents the average values of the three measures obtained from 200 random partitions of data. In most cases, the LASSO shows the best prediction performance but selects the most variables. The TLP, MCP, and SCAD show similar prediction performances and they substantially selects fewer variables than the LASSO. Among the concave penalized methods, the SCAD selects more variables than other methods. Similar to the simulation results, the methods based on regression error exhibit better prediction performances. For model size, the LASSO and SCAD based on the classification error produce more sparse models while the methods based on regression error produce more sparse model for TLP and MCP. These results suggest that the concave penalized LDA method could be a good alternative when we wish to construct the sparse model without losing the prediction accuracy much.

In this paper, we studied the high-dimensional LDA based on the concave penalized linear regression. We proved that an oracle property holds uniformly on a class of concave penalties, including the TLP, SCAD, and MCP as examples. The primary advantage of the concave penalized approach lies in its superior selection performance compared to the convex penalized approach, as supported by the simulation studies. In addition, we found that the tuning parameter selection criteria based on the prediction errors work to some extent, and hence, we may use the criteria in practice. However, we believe that there are better alternatives, such as the Bayesian information criterion, which has been proven to be useful for other penalized approaches, including the penalized linear regression (Fan and Tang, 2012).

Let ||A₂|| = sup_||u||2=1 ||Au||₂ denote the spectral norm of a matrix A. Let δ̂ = θ̂ − θ, Δ̂ = Ω̂ − Ω, and Γ̂ = Ω̂⁻¹ − Ω⁻¹.