Linear discriminant analysis (LDA) requires an estimation of the inverse conditional covariance matrix where the pooled sample covariance matrix is a manageable solution. However, the pooled sample covariance matrix is singular when the model is high-dimensional where the number of predictive variables exceeds the sample size. The performance of the LDA is undermined (Krzanowski et al., 1995) and asymptotically no better than the random guessing (Bickel and Levina, 2004) without resolving the singularity. Many works of literature have focused on this issue and suggested other alternatives by directly modifying the singular pooled sample covariance matrix or constructing relevant shrunken centroid means (Guo et al., 2006; Fan and Fan, 2008; Wu et al., 2009; Cai and Liu, 2011; Witten and Tibshirani, 2011; Clemmensen et al., 2011; Shao et al., 2011).

Meanwhile, there is a direct connection between the LDA and linear regression (Hastie et al., 2009) since the direction vector of the LDA can be obtained by using the least square estimation (LSE) (Hastie et al., 2009). However, this connection is lost when the pooled sample covariance matrix is singular. This motivated Mai et al. (2012) to study penalized LDA with the least absolute shrinkage and selection operator (LASSO) (Tibshirani, 1996). Mai et al. (2012) proved that LASSO penalized LDA finds all the features that significantly contribute to the classification when the model is high-dimensional and sparse, which is a unique theoretical work for the high-dimensional penalized LDA. Mai et al. (2012) provided various numerical studies to confirm that LASSO penalized LDA performs better than other competitors.

In this paper, we propose the penalized LDA with the moderately clipped LASSO (MCL) (Kwon et al., 2015) as an alternative to LASSO. The LASSO has been known to have higher prediction accuracy than minimax concave penalty (MCP) (Zhang, 2010) for finite samples since the shrinkage effect can increase prediction accuracy (Efron and Morris, 1975; Casella, 1985; Zhang and Huang, 2008). However, LASSO has proved to select unnecessary predictive variables even for low-dimensional linear regression (Zou, 2006). Meanwhile, MCL was designed to recover the ideal performance of LASSO by indexing a class of non-convex penalties from LASSO to MCP. Therefore MCL can select relevant predictive variables as MCP while keeping the prediction accuracy of LASSO.

We proved that, with probability tending to one, MCL penalized LDA is the same as the oracle LASSO that is a theoretically optimal LASSO obtained by using relevant predictive variables only. The equivalence can be obtained using LASSO penalized LDA as proved by Mai et al. (2012). However, MCL does not require the Strong Irrepresentable condition (Zhao and Yu, 2006) on the marginal covariance matrix, and represents a potential theoretical advantage of MCL over LASSO. We provided various numerical studies to show how MCL performs with finite samples as a better alternative to LASSO for the high-dimensional penalized LDA.

The rest of the paper is organized as follows. Section 2 introduces the penalized LDA. Section 3 introduces MCL and related statistical properties. Section 4 shows the results of numerical studies. Relevant proofs are given in the Appendix.

2.1. Linear discriminant analysis and least square estimation

Let X ∈ ℝ^p be a p-dimensional random vector of predictive variables and C ∈ {1, 2} a random class label to be identified. Given X = x, the Bayes classifier becomes

which minimizes the misclassification error, P(C ≠ φ(X)), over the set of classifiers, φ : ℝ^p → {1, 2}. The LDA (Fisher, 1936) assumes that

where μ_c and ∑ denote the mean vector and covariance matrix of the p-dimensional normal distribution. Then the Bayes classifier is equivalent to the Bayes discriminant rule: the class of C given X = x becomes 2 if x satisfies

where π_c = P(C = c), c ∈ {1, 2} are the class probabilities and

is the Bayes direction vector with θ = μ₂ − μ₁.

Let c_i ∈ {1, 2} and x_i = (x_i1, . . . , x_ip)^T ∈ ℝ^p, i ≤ n, be n samples of the class label and predictive variables. For each class label c ∈ {1, 2}, let ∑̂_c = ∑_ci=c(x_i − μ̂_c)(x_i − μ̂_c)^T /(n_c − 1) be the sample covariance matrix, μ̂_c = ∑_ci=cx_i/n_c the sample mean vector, and $n_c={\mathrm{\Sigma }}_{i=1}^{n}I(c_i=c)$ the number of samples with class label c. The LDA estimates the Bayes direction vector with the LDA direction vector:

where ∑̂ = ∑_c∈{1,2}(n_c − 1)∑̂_c/(n − 2) is the pooled sample covariance matrix and θ̂ = μ̂₂ − μ̂₁. These arguments lead to the linear discriminant rule: the class of C given X = x becomes 2 if x satisfies

where the class probabilities are estimated by the sample class proportions, π̂_c = n_c/n, c ∈ {1, 2}.

There is an intimate connection (Hastie et al., 2009) between the LDA and LSE when p ≤ n:

for some constant c > 0, where

ary and y_i = (−1)^c_in/n_ci, i ≤ n. Hence the linear discriminant rule in (2.3) is the same as:

which implies that we can cast the LDA into the framework of the LSE.

2.2. Penalized linear discriminant analysis and least square estimation

The equations in (2.2) and (2.4) fail to hold when p > n since the pooled sample covariance matrix is singular, which raises a challenging problem of estimating the Bayes direction vector. For the problem, Mai et al. (2012) proposed to estimate the Bayes direction vector by using LASSO:

for some λ > 0, where the solution can be defined even when p > n. Let Z = (Z₁, . . . , Z_p) be the centered design matrix for the LSE in (2.5), where Z_j = (I − Π₁)X_j with X_j = (x_{1 j}, . . . , x_{n j})^T, j ≤ p, and Π₁ is the projection matrix onto 1 = (1, . . . , 1)^T ∈ ℝⁿ. Then we can see that the solution can be obtained by

In addition, Mai et al. (2012) constructed an optimal discriminant rule for the penalized estimation: the class of C given X = x becomes 2 if x satisfies

whenever θ̂^Tβ̂^λ > 0. Note that the classification rule in (2.8) can be used for any linear classifier, see Proposition 2 in Mai et al. (2012) for some details. For example, the optimal discriminant rule in (2.8) reduces to the linear discriminant rule in (2.6) when p ≤ n and λ = 0 since β̂^λ = β̂^LSE = c∑̂⁻¹θ̂ implies β̂^λT∑̂ β̂^λ/θ̂^Tβ̂^λ = c.

3.1. Definition

A natural extension for LASSO in (2.7) is the non-convex penalized estimation. For example, we can use the MCP (Zhang, 2010) that has been preferred for variable selection than LASSO due to the oracle property (Fan and Li, 2001; Kim et al., 2008; Zhang, 2010). However, LASSO has been preferred for prediction than MCP since shrinkage effect often produces better prediction accuracy (Efron and Morris, 1975). As an alternative, we propose to use MCL (Kwon et al., 2015) for the penalized LDA:

where

and J_λ,γ is MCL that satisfies J_λ,γ(0) = 0 and

for some a > 1 and λ ≥ γ ≥ 0.

Note that ∇J_λ,γ(t) = λ − t/a for t < a(λ − γ) and ∇J_λ,γ(t) = γ for t ≥ a(λ − γ). Hence, MCL becomes a smooth interpolation between MCP and LASSO with two tuning parameters λ ≥ 0 and 0 ≤ γ ≤ λ. The two tuning parameters of MCL play different roles. First, λ controls the concavity of MCL near the origin for a fixed γ as it does in MCP. In the right panel of Figure 1, we can see that the concavity of MCL increases as λ does. Second, γ regularizes the amount of shrinkage in large non-zero regression coefficients for a fixed λ as it does in LASSO, which is illustrated in the left panel of Figure 1. These imply that MCL can control the sparsity and shrinkage effect simultaneously. Hence, we can expect an estimator that balances between MCP and LASSO for given finite samples as studied in Kwon et al. (2015) for high dimensional linear regression.

3.2. Asymptotic properties

In this subsection, we provide some asymptotic properties for MCL. We assume that there is a nonempty subset such that ${\beta }_j^{\text{Bayes}}\ne 0$ for and ${\beta }_j^{\text{Bayes}}=0$ for . The main results imply that MCL is asymptotically equivalent to a theoretical estimator, the oracle LASSO:

for some γ ≥ 0. Note that the oracle LASSO is simply the oracle LSE (Fan and Li, 2001) if γ = 0. The oracle LASSO is not available in practice since is unknown. However, the oracle LASSO plays an important role in developing asymptotic properties of MCL as studied in Mai et al. (2012) and Kwon et al. (2015). We can observe similar frameworks in other studies on the non-convex penalized estimation (Kim et al., 2008; Zhang, 2010; Kim and Kwon, 2012).

Before proceeding, we define some notations. For any vector a = (a_k, k ≤ s) and subset , let ${\Vert \mathbf{a}\Vert }_1={\mathrm{\Sigma }}_{k=1}^s\mid a_k\mid ,\hspace{0.17em}{\Vert \mathbf{a}\Vert }_2={({\mathrm{\Sigma }}_{k=1}^s{a}_k^2)}^{1/2}$, ||a||_∞ = max_k≤s |a_k|, supp(a) = {k : a_k ≠ 0} and . For any matrix A = (A_{k j}, k ≤ s, j ≤ t) and subsets and , let λ_min(A) be the minimum eigenvalues of symmetric A, ${\Vert \mathbf{A}\Vert }_{\infty }={\text{max}}_{k\le s}\hspace{0.17em}{\mathrm{\Sigma }}_{j=1}^t\left|A_{kj}\right|$, , and the cardinality of .

We first introduce a lemma that gives sufficient conditions for the uniqueness of a minimizer of Q_λ,γ when p ≤ n. Let Ξ_λ,γ be the set of all local minimizers of Q_λ,γ and ρ = λ_min(Z^TZ/n).

In this section, we present the results of numerical studies including simulations and real data analysis. We obtained all the penalized estimators using R package ncpen (Kim et al., 2020) that was developed based on the convex-concave procedure (Yuille and Rangarajan, 2002) and coordinate descent algorithm (Mazumder et al., 2011).

4.1. Simulation studies

The simulation studies were based on the conditional distribution in (2.1) with various scenarios. We considered two different class probabilities: π₁ = 0.5 (balanced case) and π₁ = 0.67 (unbalanced case). Given the class label, we also considered two different conditional covariance matrices: ∑ = ∑⁽¹⁾ and ∑ = ∑⁽²⁾, where ${\mathrm{\Sigma }}_{jk}^{(1)}=I(i=j)$ (no correlation) and ${\mathrm{\Sigma }}_{jk}^{(2)}={0.5}^{|j-k|}$, j, k ≤ p (power decaying correlation). For the conditional mean vectors, we set μ₁ = 0_p and μ₂ = ∑β^Bayes, where ${\beta }_j^{\text{Bayes}}=\alpha {(-1)}^{j}I(j\le q)$, j ≤ p and α was set for the Bayes misclassification error rate to be 0.2 for all scenarios.

We set n ∈ {300, 600, 900}, p ∈ {1000, 2000}, and q ∈ {5, 10} to compare finite sample performance of MCL with LASSO, including MCP as a special case of MCL. We considered three cases for MCL: MCL_k with a = 2.1 and γ = kλ̂^opt, k ∈ {1, 2, 3}, where λ̂^opt is the best tuning parameter value for LASSO (Kwon et al., 2015). We also considered three cases of MCP: MCP_k with a = k + 0.1, k ∈ {1, 2, 3} and γ = 0.

For each method, we used n training samples for direction vector estimation and n independent validation samples for selecting the tuning parameter λ. We checked the selection performance by counting the number of true positive selection (TPS), the number of false positive selection (FPS) and indicator of correct model identification (CMI): $\text{TPS}={\mathrm{\Sigma }}_{j\le q}I({\widehat{\beta }}_j\ne 0,{\beta }_j^{\text{Bayes}}\ne 0),\text{FPS}={\mathrm{\Sigma }}_{j>q}I({\widehat{\beta }}_j\ne 0,{\beta }_j^{\text{Bayes}}=0)$, and CMI = I(TPS = q, FPS = p − q). Furthermore, we compared the misclassification error rate (ERR) obtained from 2n independent test samples by using the classification rule in (2.8). We repeated each simulation 200 times and summarized the averages of the four measures in Table 1 and Table 2, including graphical illustrations in Figure 2 and Figure 3.

Table 1 shows the results when ∑ = ∑⁽¹⁾. First, TPS converged to q as n increased for all methods while LASSO performed best for each fixed n. Second, FPS decreased as n increased only for MCL and MCP, leading to an increase in CMI. Hence, our theoretical results were supported by the results since MCL and MCP correctly fitted the model as n increased while LASSO overfitted the model regardless of sample size. These trends continued even under the class imbalance and dimension increase. Third, ERR tended to decrease as n increased for all methods. Especially MCL₁ attained the lowest ERR for each fixed n and was followed by MCL₂ and LASSO.

Table 2 shows the results when ∑ = ∑⁽²⁾, where MCPs performed the best for all measures. All selection performance measures tended to be undermined but the patterns were similar in that MCL and MCP fitted the model correctly while LASSO overfitted. For ERR, MCPs performed the best and MCL₁ ranked the second for all cases.

We conclude that MCL can be a nice alternative to LASSO for the high-dimensional penalized LDA. The MCL can correctly identify the sparse Bayes direction vector, keeping almost the same prediction accuracy as LASSO. The MCL₁ performed quite well regardless of the simulation designs considered in this paper, which aligns with the recommendation for the heuristic choice of γ = λ̂^opt in the linear regression (Kwon et al., 2015).

4.2. Analysis of micro-array samples

The R package datamicroarray (John, 2016) provides a collection of high-dimensional microarray data sets. We chose four data sets (Burczynski et al., 2006; Chin et al., 2006; Chowdary et al., 2006; Gordon et al., 2002) to illustrate how MCL performs with real samples. We first applied the sure independence screening procedure (Fan and Song, 2010) so that the first top d ∈ {400, 800, 1600} predictive variables with the largest marginal regression coefficients were used for the penalized LDA. For comparison, we applied the leave-one-out cross-validation procedure for each data set and calculated the number of incorrectly classified samples (error) as well as the average of the numbers of non-zero regression coefficients (sizes), where the tuning parameters were selected by the 10-fold cross-validation procedure.

Table 3 summarizes the results. In most cases, LASSO had the best prediction accuracy selecting the most variables. The MCPs had the worst prediction accuracy but selected the least variables. The best prediction accuracy occurred when d = 400 for all cases where MCL₁ had the same prediction accuracy as LASSO while selecting fewer variables than LASSO. In addition, MCL₂ selected fewer variables than LASSO without losing prediction accuracy much.

In this paper, we studied the high dimensional penalized LDA with MCL. The nature of MCL produces the same shrinkage effect as LASSO and makes the selection process the same as MCP. Therefore MCL shows similar or better prediction accuracy than LASSO correctly recovering the sparsity of the direction vector. We proved that MCL is selection consistent under reasonable regularity conditions, which was supported by various numerical experiments. One disadvantage of MCL compared with LASSO may be the additional tuning parameter γ. But the heuristic choice of γ = γ̂^opt or γ = 2γ̂^opt performed well as shown in the numerical studies. Further research on MCL could focus on the theoretical justification of the choice of γ, which is not addressed in this paper.