LDA via L1-PCA of Whitened Data

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 225 LDA via L1-PCA of Whitened Data Rubén Martín-Clemente , Member, IEEE, and Vicente Zarzoso , Senior Member, IEEE Abstract—Principal component analysis (PCA) and Fisher’s linear discriminant analysis (LDA) are widespread techniques in data analysis and pattern recognition. Recently, the L1-norm has been proposed as an alternative criterion to classical L2-norm in PCA, drawing considerable research interest on account of its increased robustness to outliers. The present work proves that, combined with a whitening preprocessing step, L1-PCA can perform LDA in an unsupervised manner, i.e., sparing the need for labelled data. Rigorous proof is given in the case of data drawn from a mixture of Gaussians. A number of numerical experiments on synthetic as well as real data confirm the theoretical findings. Index Terms—Fisher’s linear discriminant analysis, L1-norm, principal component analysis. I. INTRODUCTION FISHER’S linear discriminant analysis (LDA) and principal component analysis (PCA) can be considered as two pillars ofdataanalysis[1].Givenadatasetwithtwoclasses,LDAfindsa projectionof thedata pointsonto aone-dimensional spacewhere theclassesarewellseparated.WhenthetwoclassesareGaussian with equal covariances, LDA yields the optimal Bayes classifier. LDA is a supervised technique, whose performance depends heavily on the availability of correctly labelled data [2]. On the other hand, PCA finds the linear projection best fitting the data in the least-squares sense, which is also the projection of maximum variance of the dataset. As a fully data-driven technique, PCA is unsupervised and does not require labelled samples [3]. Spurred by the versatility of PCA, the last decade has witnessed a flurry of research into alternative criteria aimed at enhancing its capabilities or alleviating its limitations in various operating conditions. One such criterion, in particular, is based on replacing the L2-norm of classical PCA by the L1-norm, which is more robust against outliers, thus giving rise to L1- PCA [4]. This technique, compared to other robust versions of PCA [5]–[8], is particularly intuitive and simple, as well as invariant to the rotation of the data, all of which justifies its growing interest and use in a wide range of applications, e.g., sensor-array processing, image fusion, video surveillance, robust face recognition, et cetera [9]–[13]. In addition, a number Manuscript received October 13, 2018; revised July 14, 2019; accepted November 13, 2019. Date of publication November 25, 2019; date of current version December 27, 2019. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Yuantao Gu. This work was supported by the Spanish Ministry of Economy and Competitiveness under project TEC2017-82807-P. The work of R. Martín-Clemente was supported by professorship at the I3S Laboratory in July 2018 funded by an IFA scholarship fromUniversitéCôted’Azur. (Corresponding author: Ruben Martin-Clemente.) R. Martín-Clemente is with the Departamento de Teoría de la Señal y Comunicaciones, Escuela Superior de Ingeniería, Universidad de Sevilla, Seville 41092, Spain (e-mail: [email protected]). V. Zarzoso is with the Université Côte d’Azur, CNRS, I3S Laboratory, Sophia Antipolis Cedex 06903, France (e-mail: [email protected]). Digital Object Identifier 10.1109/TSP.2019.2955860 of efficient L1-PCA algorithms have been proposed in the recent literature (see e.g. [14]–[18]). Even though L1-PCA and LDA are apparently disparate techniques with very different purposes, the present contribution shows that a link actually exists between both approaches. We prove that, for whitened (or sphered) data, L1-PCA can perform LDA in an unsupervised fashion, that is, sparing the need for training data. This result is of theoretical interest and opens interesting research perspectives for performing LDA using L1- PCA algorithms. A number of numerical experiments validate the theoretical findings in a variety of simulation scenarios. The paper is organized as follows. Section II reviews the basics of Fisher’s LDA and PCA and other related approaches such as kurtosis optimization. Section III analyzes the L1-LDA criterion using a Gaussian assumption, and shows it to perform LDA in an unsupervised manner under certain prescribed conditions. Links with other techniques are also established and working algorithms reviewed. Generalization to more than two classes is addressed in Section IV. Computer experiments validating the theoretical findings on synthetic and real data are reported in Section V. The concluding remarks of Section VI bring the paper to an end. The notation used in the paper is as follows. Lightface letters (e.g., a,A) represent scalar quantities, which may be functions, univariate random variables, deterministic constants or indices. The distinction will be clear from the context. Boldface lowercase (x) and uppercase (X) letters, respectively, stand for vectors and matrices, either deterministic or random depending on the context. II. BACKGROUND This section reviews the basic concepts on LDA and PCA, recalling some connections with other techniques. This material will be useful for the analysis in Section III. A. Assumptions The following conditions are assumed throughout the paper. Consider p-dimensional observations x∈Rp, where each observation xbelongs to one of two classes C1and C2. Suppose that the two classes can be described by probability density functionsf1andf2,respectively.ThepriorprobabilityofclassCi is denoted by πi,i=1,2,with π1+π2=1. Distribution fihas a mean µi=E{x|x∈C i}and an invertible positive definite covariance matrix Vi=Cov(x|x∈C i)=E{(x−µi)(x− µi)|x∈C i},i=1,2, where E{·}is the expectation operator and (·)stands for transpose. The following typical assumption can be made without loss of generality: Assumption A1: The data have zero mean, i.e., µ=E{x}= 2 i=1 πiµi=0. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 226 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 Under this assumption, one can easily show that: µ1=−π2Δµ(1a) µ2=π1Δµ(1b) where Δµ=(µ2−µ1). B. LDA Formulation LDA [1] aims at projecting the data along a direction a∈Rp, resulting in the projection y=ax.(2) A reasonable criterion is to select a vector asuch that the corresponding projections yobtained from C1and C2are separated as much as possible. To this end, the derivation of LDA involves the so-called within-class scatter matrix, defined as SW= 2  i=1 πiVi(3) which quantifies the average variation of the data in the classes. Similarly, the between-class scatter matrix SB= 2  i=1 πiµiµ i=π1π2ΔµΔµ(4) contains the distance between the class centers. Equation (1) has been invoked to reach the last term of Eqn. (4). To find direction a, LDA maximizes the Rayleigh quotient: JLDA(a)= aSBa aSWa(5) which is simply obtained at [1] SWa∗=Δµ.(6) The rationale behind LDA cost (5) can be explained by the following: observe that the mean and variance of the projected data of class Cican be expressed, respectively, as mi=aµiand σ2 i=aViai=1,2.(7) Then simple algebra shows that JLDA(a)= π1π2(Δm)2 π1σ2 1+π2σ2 2 with Δmdef =m2−m1=aµ2−aµ1=aΔµ.(8) From here it is apparent that LDA searches for a vector a such that: (i) m1and m2are far apart from each other, and (ii) the spread of the projections around m1and m2is small. Both conditions try to prevent the overlapping and maximize the separation of the projected classes, thereby easing their discrimination. To compute class statistics {µi,Vi}2 i=1, LDA needs labelled data — also known as training data — where the class to which each observed vector xbelongs is known. This is why LDA is a supervised classification technique. 1) Whitened Data: LDA admits an interesting interpretation when the classes are whitened or sphered, as explained next. Let us now consider the case in which the classes are homoscedastic, i.e., they have the same covariance matrix, V1=V2=V, which is a commonly accepted assumption in LDA. Then, Eqn. (3) reduces to SW=Vand Eqn. (6) simplifies into Va ∗=Δµ.(9) The projection onto vector (9) gives y=xa∗=ΔµV−1x. It is now useful to calculate the eigendecomposition of the within-class covariance matrix, given by V=UDU. Then we have V−1=UD−1U, and y=Δ ˜ µ˜ x where ˜ x=D−1/2Uxand Δ˜ µ=D−1/2UΔµ. This transformation basically whitens (or spheres) the data because the covariance of ˜ xis the identity matrix for both classes, Cov(˜ x|x∈ Ci)=D−1/2UVUD −1/2=I,i=1,2. Therefore, LDA is equivalent to whitening each class and then projecting the data along the line joining the whitened class centroids, defined by vector Δ˜ µ. 2) Kurtosis-Based Unsupervised LDA: The kurtosis, defined as the fourth-order central moment divided by the squared variance, is a measure of the gaussianity, peakedness and bimodality of a distribution [19]. We recall it here because the following theorem, reproduced from [20] using our notation, shows that LDA is also closely related to this statistic. Theorem 1: Let xbe a p-dimensional random variable distributed as π1f1(x)+π2f2(x), where π1+π2=1 and fi, i=1, 2, is a Gaussian distribution with mean µiand covariance matrix Vi=V, the same for both distributions. Let abe a unit-norm vector on Rpand y=ax.Ifasatisfies Va=Δµ(10) then it maximizes the absolute kurtosis of y. Furthermore, these directions minimize the kurtosis if |π1−1/2|<1/√12, and maximize it otherwise. A related result for the case of different covariance matrices can be also found in [20]. It is interesting to observe that (10) is equivalent to (9) but can be obtained without prior knowledge of the data allocation, i.e., without the need for training data. This result opens the way for kurtosis-based unsupervised classification techniques [21]. The main drawback of this approach is that kurtosis is very sensitive to outliers since, by raising the projections to the fourth power, the effects of the data points far from the nominal distribution can easily overshoot, leading to poor results in the presence of faulty data. 3) Generalized LDA: Finally, the LDA solution (6) can be also generalized to weighted within-class scatter matrices of the form: ¯ SW= 2  i=1 πiβiVi(11) where βi,i=1, 2, are constants such that ¯ SWis positive or negative definite. Such weighted scatter matrices lead to admissible LDA-type classifiers for multivariate normal distributions MARTÍN-CLEMENTE AND ZARZOSO: LDA via L1-PCA OF WHITENED DATA 227 with different covariance matrices [22]. Notice that the L1-norm based criterion studied in Section III will be shown to be a particular generalized LDA classifier under certain conditions. C. L2-PCA and L1-PCA In its original formulation [3], classical PCA seeks directions maximizing the variance of the projection (2) according to the quadratic criterion: max a E{y2}subject to a2=1.(12) This technique is also known as L2-PCA because, when operating on finite sample observations, the expectation in (12) is usually replaced by the L2-norm of the vector of projected samples.PCAwasoriginallydevelopedfordimensionalityreduction of multivariate data, in general. The solution to problem (12) can simply be computed as the dominant eigenvector of the overall data covariance matrix, which under the assumptions of Section II-A is given by V0=E{xx}=SW+SB.(13) If the classes are well separated, the between-class scatter typically dominates the within-class scatter, i.e., SWFro SBFro (14) where ·Fro represents the Frobenius norm, and then: V0≈SB=π1π2ΔµΔµ.(15) As a result, L2-PCA defines a projection in the direction joining the class centroids, Δµ. This result reminds us of the LDA solutionfor whitenedclassesas presentedin theprevioussection but, in contrast to LDA, only in certain situations it provides a good direction of projection for separating the two classes. For the interested reader, a comparison between LDA and PCA in image recognition tasks is given for example in [23]. Though mathematically appealing, the L2-norm is rather sensitivetoimpulsivenoiseoroutlierssincesquaringtheprojections overemphasizes the effects of far-off data points. As presented in [4], a natural extension of L2-PCA to deal with this drawback is given by max a E{|y|}subject to a2=1.(16) The absolute value in Eqn. (16) is replaced in practice by the L1-norm of the vector of projected samples and thus the above criterion is referred to as L1-PCA. Note that the above two PCA criteria are not equivalent in general, though they can provide similar results when the data contain no outliers [4], [16]. In addition, note also that none of them require labeled training data. III. LINK BETWEEN L1-PCA AND LDA Inspired by the results reviewed in the previous section, we consider the following unsupervised criterion for linear discriminant analysis: max a E{|y|}subject to E{y2}=1 (17) with y=axand hence E{y2}=E{(ax)2}=aV0a.Observe that the constraint differs from that used in L1-PCA (a2=1) but still prevents E{|y|}from increasing by a mere increase of the magnitude of a. In the rest of the paper, the above criterion is referred to as L1-norm based unsupervised LDA (L1-uLDA). The L1-uLDA criterion is analyzed in the remainder ofthissection.Inparticular,theunit-varianceconstraintisshown to be the ingredient confering discriminative capabilities. The L1-uLDA criterion is quite general in that it does not require the data to be prewhitened. Without limiting the foregoing, it is well-known that the unit-variance constraint equates to a unitnorm constraint for whitened data, so that L1-uLDA becomes L1-PCA in this case. This observation enables us to make the link between L1-PCA and LDA explicit in Section III-D. It is concluded that, for whitened data, L1-PCA can perform LDA in an unsupervised fashion, without information about class labels. For unwhitened data, which is the more general case, this result still holds since a simple sphering or whitening preprocessing step allows us to carry out L1-uLDA by L1-PCA. A. Justification of the Criterion Anintuitiveexplanationof how thecriterion(17) worksis that the absolute value can be considered as a measure of distance from the origin; as a result, it seeks a direction maximizing the average distance from zero of the projected data points, which in the end, thanks to the zero-mean assumption, favours the emergence of two distinct clusters on either side of the origin. The second power used in L2-PCA (12) can also be considered as another measure of distance from the origin but it leads to the direction of maximum variability of the data, which, as it is well-known, does not necessarily produce linear discriminant results. The main point is that L1-PCA is endowed with discriminative power thanks to the unit-variance constraint. We will show that, after whitening the data, L1-PCA is able to find the direction of the line joining the class centroids. Furthermore, as previously shown in Section II-B2, there exists a close connection between LDA and the optimization of the kurtosis of the linearly projected data. An additional motivation for criterion (17) lies in the fact that, as recently shown in [24], the kurtosis is related to the mean of the absolute value of the random variable. To see this, let us standardize the random variable yto have zero mean and unit variance and, assuming that moments exist up to order four, consider the fourth-order Gram-Charlier expansion of its probability density function [25]: f(y)≈exp(−y2/2) √2π1+λ3 6(y3−3y)+κ4 24(y4−6y2+3)  where λ3=E{y3}and κ4=(E{y4}−3) denotes the excess kurtosis (i.e., the kurtosis minus three) of y. Then, using that ∞ 0yuexp(−y2/2)dy =0.5−((u−1)/2)Γ((u+1)/2),Γ(·)being the Gamma function, some algebra shows that E{|y|} =∞ −∞ |y|f(y)dy=2 π1−κ4 24.(18) As a result, maximizing (resp. minimizing) E{|y|}is equivalent to minimizing (resp. maximizing) the kurtosis. Although one may question whether the Gram-Charlier expansion can represent adequately the distribution of the data, this equivalence 228 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 Fig. 1. Plot of function g(α)defined by eqn. (23). The constant 1/√πand the straight line y=αare represented by the dotted and dashed-dotted lines, respectively. together with Theorem 1 (Section II-B2) suggest exploring the link between the averaged absolute value and LDA. B. Theoretical Analysis in the Gaussian Case As it is usual in the study of LDA, we assume that the clusters arenormallydistributed,i.e.,x|Ci∼N(µi,Vi),i=1,2.This is a probabilistic model that has been extensively applied in signal processing. It follows that the observations are distributed as a mixture of Gaussians, x∼2 i=1 πiN(µi,Vi). Hence, any linear projection as in Eqn. (2) will also be distributed according to a Gaussian mixture model y∼2 i=1 πiN(mi,σ2 i), with probability density function f(y)= 2  i=1 πi 2πσ2 i exp −(y−mi)2 2σ2 i(19) where miand σ2 iare defined in (7). As a consequence of Assumption A1, the mean value of yis also zero. The objective function and the constraint can be expressed as follows. The second-order moment of yis readily given by E{y2}=aV0a(20) where V0is the covariance matrix of xdefined in (13). Next, using the change of variable x=(y−mi)/σiand letting αi=mi √2σi (21) some algebraic manipulations show that E{|y|} =∞ −∞ |y|f(y)dy=√2 2  i=1 πiσig(αi)(22) with g(αi)= 1 √πexp(−α2 i)+αierf (αi)(23) where the error function is defined as erf(z)= 2 √πz 0e−t2dt.It holds that (23) is an even function and that lim|α|→∞ g(α)=|α|. Fig. 1 plots this function. Appendix A shows that the stationary point of the constrained optimization problem (17) is given by ¯ Sa =δΔµ(24) where ¯ S= 2  i=1 πiβiVi(25a) βi=2/π σi exp(−α2 i)−E{|y|} i=1,2(25b) δ=π1π2(γ2−γ1)(25c) γi=miE{|y|}−erf(αi)i=1,2(25d) and αi,i=1,2, were defined in Eqn. (21). Furthermore, it is important to point out that solution (24) locally maximizes criterion (17) when |αi|,i=1, 2, are sufficiently large (see Appendix A). To shed some light on this result, some comments are in order. Firstly, observe that when |α1|and |α2|are large, the magnitude of the projected means |mi|are large relative to the standard deviations σi. Therefore, the projected data form clusters and, because of the zero-mean assumption, one cluster is to the left and the other to the right of the origin, which is the desirable outcome in linear discrimination [22]. Secondly, when |αi|is sufficiently large, exp(−α2 i)≈ 0and erf(αi)≈sign(αi). In that case, and recalling that lim|α|→∞ g(α)=|α|, expression (22) simplifies to E{|y|} ≈ π2|m2|+π1|m1|=π2m2−π1m1 =2π1π2aΔµ=2π1π2Δm(26) where we have assumed, without loss of generality, that m1< 0<m 2. Observe that maximizing (26) with a unit-norm constraint, as in the original L1-PCA criterion (16), would yield a projection on the direction of the line joining the class centroids, Δµ, which is the direction of maximum variance: Section II-C recalled that, under condition (14), L2-PCA also finds the direction Δµ. It is the unit-variance constraint in L1-uLDA variant (17) that endows L1-PCA with LDA’s discriminative power, as shown throughout this paper. C. Link With LDA Now we are ready to establish the formal connection with LDA. In general, if β1and β2in Eqns. (24)–(25) are both positive (resp. negative) then ¯ Sis positive (resp. negative) definite because V1and V2are both positive definite. In this case, the L1-uLDA solution (24)–(25) is a particular instance of LDA classifier (6) with generalized covariance (11). The only additional requirement for (24) to define an LDA solution is that δ=0, which, according to Eqns. (25c)–(25d), is equivalent to the constraint: Δm=erf(α2)−erf(α1) E{|y|} .(27) Next, suppose that the projected data form clusters where the meanis largerelative tothe standarddeviation,in sucha waythat |αi|,i=1,2, are sufficiently large (asymptotic regime). Hence, the exponentials exp(−α2)vanish in (25b), yielding β1=β2≈ MARTÍN-CLEMENTE AND ZARZOSO: LDA via L1-PCA OF WHITENED DATA 229 −E{|y|}. Then, the generalized covariance matrix (25a) can readily be expressed as ¯ S≈−E{|y|} 2  i=1 πiVi∝SW which leads to Fisher’s LDA solution (6). Taking into account Eqn. (26), the condition on δnow becomes: Δm=1 √π1π2 (28) where we have supposed that m1<0<m 2without loss of generality. Appendix B proves that condition (28) is always satisfied. Thediscussionsimplifieswhenthe classesarehomoscedastic, i.e., V1=V2=V. Then, the L1-PCA scatter matrix (25a) reads: ¯ S=2  k=1 πkβkV∝V so that Eqn. (24) always turns out to be equivalent to Eqn. (9) defining LDA in the homoscedastic scenario as recalled in Section II-B. This result also coincides with Theorem 1 in Section II-B2 for the kurtosis-based unsupervised classifier, and provides additional evidence of the relation between the proposed unsupervised criterion and LDA. The above analysis can be summarized as follows: i) The critical points of L1-uLDA criterion (17) present the general form of solutions (6) with general covariance given by Eqn. (25). ii) The equivalence between L1-uLDA and LDA is exact if the projection onto the directions they define produces two sufficiently separated clusters on both sides of the origin so that |αi|,i=1,2, are sufficiently large. iii) Under this condition, LDA is obtained in an unsupervised fashion by maximizing the proposed L1-uLDA criterion. Note that the condition guaranteeing the equivalence between LDA and L1-uLDA is mild, as it coincides with the condition for LDA to yield proper discrimination results. A complementary discussion on the conditions under which |αi|,i=1,2,can be considered sufficiently large is given in Appendix C. Despite the theoretical equivalence in ideal conditions, L1-uLDA operates in an unsupervised fashion. This theoretical analysis will be put to the test in the experimental analysis of Section V. D. Algorithm 1) Derivation: We haveseenthattheunit-varianceconstraint in (17) is crucial for L1-uLDA to yield discriminative solutions. However, classical algorithms for L1-PCA work under a unitnorm constraint, as in Eqn. (16). A sphering or whitening preprocessing step transforms the variance constraint into the norm constraint, thus enabling the use of classical algorithms for L1- PCA. Given V0, the covariance matrix of x, and its eigenvalue decomposition (EVD) V0=QD0Q, consider the change of variable z=D−1/2 0Qx. This transformation spheres or whitens the data, as the covariance matrix of zbecomes the identity.Let w=D1/2 0Qaandobservethat y=wz=ax, implying that E{y2}=wE{zz}w=ww=w2 2. As a result, for whitened data the proposed criterion (17) can be rewritten as max w E{|y|}subject to w2=1 (29) where y=wz. Hence, the unit-power constraint on the projected data is transformed into a unit-norm constraint on the projection vector. The constrained optimization problem (29) is the basis of L1-PCA and has been studied in a number of recent works [4], [16]. In particular, it can be shown that the fixed point iteration 1) y=w nz 2) wn+1 =E{zsign(y)} E{zsign(y)}2 monotonicallyincreasestheabsolutevaluecriterionE{|y|}after each iteration [4], and so the algorithm converges at least to a local maximum — similarly, given a set of unlabeled data points x1,...,xTsampled from the random variable x,step2 transforms into wn+1 =T t=1ztsign(w nzt) T t=1ztsign(w nzt))2 . The most notable feature of this simple algorithm is that no parameters need to be tuned. Alternatively, one can use the polynomial time approaches in [26]–[28]. L1-PCA algorithms with guaranteed convergence to a global maximum have also been proposed, although they come at the expense of increased computational complexity [16], [29]. In addition, a simplified and faster, yet suboptimal, version of [16] can be found in [17]. Insummary,the algorithm(assumingunwhitened data)would involve the following steps: 1) Compute and subtract the mean from the data. 2) Whitening: z←D−1/2 0Qx, where V0=QD0Qis the EVD of the covariance matrix of x. 3) Solve (29) to obtain w∗. This step is the standard L1-PCA of the sphered data and any of the above mentioned algorithms (e.g., [4], [16], [17]) can be used for its computation. 4) Undo whitening: return a∗=QD−1/2 0w∗. Interestingly, this algorithm has recently been shown to also carry out independent component analysis (ICA) with increased protection against outliers [24]. In other words, if the data follow the ICA model, this algorithm can allow its estimation. Otherwise, if the data have two classes, this same algorithm can perform LDA. 2) Further Comments: After sphering, the data covariance matrixequalstheidentity,i.e.,V0=E{xx}=I.Since V0= SW+SB, it follows that SW=V0−SB=I−π1π2ΔµΔµ where we have invoked Eqns. (4) and (13). This matrix can be easily inverted using the well-known Sherman-Morrison formula, yielding: S−1 W=I−π1π2ΔµΔµ 1+π1π2Δµ2. 230 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 It follows that the optimal LDA projection is in the direction of the line joining the class means: a∗=S−1 WΔµ=Δµ 1+π1π2Δµ2. If the classes are well-separated, Δµalso coincides with the direction of maximum variance (as recalled in Section II-C). Consequently,we may thinkthat Δµ(andhencea) canbe easily estimated in an unsupervised manner by L2-PCA. However, after data sphering all directions yield the same variance and henceΔµcannot be identified from a variance criterion. In other words, L2-PCA becomes inoperative. To fix this problem, the proposed L1-uLDA algorithm replaces L2-PCA with L1-PCA. Having replaced the square by the absolute value in the objective function, L1-PCA is able to identify Δµas shown throughout this section. IV. THE MULTICLASS CASE In the preceding section of the paper, the link between LDA and the constrained L1-PCA giving rise to L1-uLDA has been established in the binary (two-class) case. This section extends theseresultstothegeneralscenariowheredatamaybecomposed of more than two classes. We begin the derivations by recalling LDA in the multiclass case. 1) Multiclass LDA: In the multiclass case, where xis drawn from one of cclasses C1,...,Cc,thebetween-class and withinclass scatter matrices are defined respectively as [30]: Sc B= c  i=1 πi(µi−µ)(µi−µ)= c  i=1 πiµiµ i Sc W= c  i=1 πiVi. where we have assumed again that the total mean vector of x is zero, i.e., µ=E{x}=c i=1πiµi=0, and maintained the previous notation: µi=E{x|x∈C i}is the mean of Ci,Vi= Cov(x|x∈C i)=E{(x−µi)(x−µi)|x∈C i}istheclass covariance matrix, and πi=P(x∈C i)is the corresponding a priori probability. Regardless of the value of c, it always holds that the data covariance matrix V0=E{xx}can be decomposed as V0=Sc B+Sc W.(30) Multiclass LDA addresses the more general problem of projecting the data onto a q-dimensional subspace, with q<p, while also retaining as much of the class discriminating information as possible for the subsequent classification stage. The LDA projection is given by y=Ax, where Ais the (p×q)matrix whose columns are the dominant eigenvectors of the matrix pencil (Sc W,Sc B), given by the solutions of the generalized eigenvalue (GEVD) problem [30]: Sc Wa=λSc Ba.(31) This can be seen as the solution of two coupled optimization problems involving quadratic forms subject to the same unit variance constraint. In the first place, because of relation (30), the GEVD problem (31) can readily be expressed as: Sc Ba=λV0a(32) where λdef =1/(1 + λ). Some algebraic manipulations prove that Eqn. (32) is the solution of the constrained problem: max aSc Basubject to E{y2}=1 (33) since E{y2}=E{(ax)2}=aV0a. Likewise, problem (31) can also be expressed as: Sc Wa=λV0a(34) with λ def =λ/(1 + λ), which is the solution of min aSc Wasubject to E{y2}=1.(35) Problems (32)–(33) and (34)–(35) are indeed inextricably intertwined because, under the unit-variance constraint: arg max a aSc Ba=argmax a a(V0−Sc W)a =argmax a (1 −aSc Wa) =argmin a aSc Wa and λ+λ =1. Therefore, the dominant eigenpairs of (32) and the maxima of (33) are directly linked to the minor (least significant) eigenpairs of (34) and the minima of (35), respectively. Finally, Eqn. (32) can be also expressed as λV0a=Sc Ba=c  i=1 πiµiµ ia= c  i=1 (πimi)µi,(36) wheremi=µ ia.Hence,we see thatthecolumns of theoptimal matrix Alie in the space spanned by µ1,...,µc. In addition, pre-multiplying Eqn. (36) by aand imposing the unit-variance constraint shows that the dominant generalized eigenvector of (32) is given by λ max =maxaSc Ba=max c  i=1 πi(mi)2.(37) Inotherwords,theLDAsolutionliesinthedirectionmaximizing the average squared means of the projected classes under the unit-variance constraint. In a totally similar fashion, it is simple to observe that the minor generalized eigenvector of (34) under the unit-variance constraint is given by λ min =minaSc Wa=min c  i=1 πiσ2 i which minimizes the average variance of the projected classes. In summary, multiclass LDA tries to find projections maximizing the average squared mean while minimizing the average variance of the projected classes under a unit-norm constraint on the projected data. Inspired by this result, we show next that L1-uLDA follows a related approach in the multiclass scenario. 2) Multiclass L1-Ulda in Low Dispersion Regime: A slight variant of LDA is obtained by replacing the squares by absolute values in Eqn. (37), yielding max c  i=1 πi|mi|subject to E{y2}=1.(38) MARTÍN-CLEMENTE AND ZARZOSO: LDA via L1-PCA OF WHITENED DATA 231 Some algebra shows that the solution fulfils V0a∝ c  i=1 πisign(mi)µi(39) which, as in classical LDA, lies in the subspace spanned by the class means [cf. Eqn. (36)]. This alternative criterion remains supervised, since computing mirequires prior knowledge of the class labels. We now prove that the L1-uLDA criterion optimizes (38) in an unsupervised fashion when the data classes are tightly concentrated aroundtheir means.Tocompletethe linkwith LDA,thenext section will analyze the other extreme case where the within-class spread is large. Recall that the L1-uLDA criterion is given by max E{|y|} subject to E{y2}=1.(40) We first consider the case where the classes are tightly concentrated around their means, i.e., trace(Vi)≈0and σi→0, ∀i=1,2,...,c. Hence, the probability density function (pdf) of the projected data accepts the approximation: fy(u)−→ σi→0 c  i=1 πiδ(u−mi). Therefore: E{|y|} =+∞ −∞ |u|fy(u)du −→ σi→0 c  i=1 πi+∞ −∞ |u|δ(u−mi)du = c  i=1 πi|mi|. The L1-uLDA criterion (40) indeed asymptotically optimizes (38) unsupervisedly, i.e., without requiring knowledge of the class labels. Remark that this result is independent of the underlying data distribution. The only requirement is that each class be closely concentrated around its mean. 3) Multiclass L1-uLDA in High Dispersion Regime: We now turn to the case where the within-class dispersion is large as compared with the between-class dispersion. To address this more involved scenario, the assumption is now made that the data xfrom class Ciare Gaussian distributed. The pdf of xis therefore x∼c i=1 πiN(µi,Vi). Under this model it can be shown that E{|y|} =∞ −∞ |y|f(y)dy=√2 c  i=1 πiσig(αi)(41) where αiand g(αi)were respectively defined in Eqns. (21) and(23),andσ2 iisthevarianceof theprojectedithclass.Aseries of algebraic manipulations, similar to those for the two-class case, yield that the stationary points of the L1-uLDA criterion verify the equation V0a=1 E{|y|} c  i=1 πi2/π σi e−α2 iVia+erf(αi)µi. (42) If the spread of the clusters is large in relation to their respective means, αiare close to zero, implying exp(−α2 i)≈1and erf(αi)≈0. Then (42) can be rewritten as the GEVD problem ¯ Sc Wa=λV0a(43) where ¯ Sc W= c  i=1 πiβiVi,with βi=2/π σi . Equation (43) is the multiclass LDA solution with a generalized weighted within-class scatter matrix. Furthermore, it becomes equivalent to the traditional LDA solution (34) when all classes have equal covariance matrices, Vi=V,∀i=1,2,...,c.On the other hand, remark that in the low dispersion regime αi become large and then exp(−α2 i)≈0and erf(αi)≈sign(αi), leading to the same conclusions as in the previous section. We have established the link between LDA and L1-uLDA in extreme cases of within-class dispersion. In the general case, the actual solutions of the L1-uLDA criterion (40) will lie somewhere between these two extremes, preserving in any case theabilityof multiclassLDAtodistinguishbetween thedifferent classes, as demonstrated by the numerical experiments reported next. V. NUMERICAL EXPERIMENTS A number of computer simulations in a variety of experimental conditions are performed to validate the theoretical study developed in this paper and to test the equivalence to LDA while sparing the need for labelled data (unsupervised operation). Other techniques for unsupervised and semi-supervised LDA have been proposed in the literature (see e.g. [31]–[33]). It should be kept in mind, however, that to limit the scope of the paper our focus is only on the connection between L1- PCA and LDA, without claiming its potential superiority over existing alternative techniques for classification or clustering. A comparative performance analysis should be the topic of future research. Note also that other complementary analyses of L1-PCA(robustness againstoutliers, lineardimensionreduction capabilities, comparison with other robust PCA approaches, etc.) can already be found in the cited literature (see e.g. [4, 17]). Our experiments include synthetically generated as well as realdata.Forsimplicity,weapplyourteststoL1-uLDA,whichis carriedoutbyapplyingL1-PCAtoprewhitened data.Toperform L1-PCA in step 3 of the L1-uLDA optimization algorithm in Section III-D, and unless otherwise stated, we employ the bit flipping algorithm presented in [17]. A free MATLAB implementation of this algorithm is provided in [34]. A. Synthetic Data Bivariate data can be easily visualized as points in a twodimensionalscatterplot.Forthisreason,weusethemtoillustrate the performance of L1-uLDA first before moving on to the general case of more than two dimensions. 1) Finite Sample Size: The theoretical analysis of the previous section relies on an ensemble or distributional characterization with the implicit assumption of infinite sample size. The first experiment evaluates the ancillary L1-uLDA criterion in short data records, and shows the existence of spurious 232 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 Fig. 2. Equivalent L1-uLDA criterion (44) computed from a realization of T= 200 points (red line) and T=50points (blue line) as a function of the angle formed by the projection direction with the horizontal axis. Spurious local maxima are clearly visible in the latter case. local maxima for insufficient sample size. We generate twodimensional random samples from a mixture of two equiprobable (π1=π2=0.5) Gaussian distributions with means µ1= −[1.5,0]and µ2=[1.5,0], and identical covariance matrix V1=V2=I, where Iis the identity matrix. Let y=a(θ)x be the one-dimensional projection of the data vector xin the direction a(θ)=[cos(θ),sin(θ)]. It can be readily shown that the L1-uLDA criterion (17) is also equivalent to maximizing the unconstrained function J(θ)=E{|y|} σy (44) where σyis the standard deviation of y. Indeed, the optimization of a Rayleigh quotient as that of (44) is equivalent to optimizing the numerator while constraining the denominator. The maxima and the minima of (17) can be easily studied by plotting J(θ). The red line in Fig. 2 plots function J(θ)estimated from a realization of T= 200 data samples. The curve passes through a global maximum at θ=0degrees, which is consistent with the fact that the optimal direction resulting from Fisher’s discriminant is parallel to the x1-axis as is easy to check (line joining the means; see Eqn. (9)). Next, the experiment is repeated, though this second time the objective function is estimated from only the first 50 samples of the previous 200-samples dataset. The corresponding objective function is shown in the blue line of Fig. 2. Observe that several spurious local maxima appear when reducing the sample size. We conclude that it is convenient to carry out L1-PCA using globally convergent algorithms such as that of [16] if the sample size is small. The bit flipping algorithm of [17], though suboptimal, achieves good convergence in all our experiments. 2) Comparison with LDA and L2-PCA: Since L1-PCA is related with L2-PCA (Section II-C), one may naturally wonder about their comparative behavior. To this end, Fig. 3 shows the scatter plot of the T= 200 data points used in the previous experiment, where data from class 1 and class 2 appear as red crossesandbluecircles,respectively.The redlinemarked‘L1’ is the direction of the optimal L1-uLDA projection. Projecting the whole dataset onto this line gives the histogram plotted in Fig. 4, clearly showing the clustered structure of the data. Similarly, the blacklinelabelledas‘PCA’inFig.3 isthedirection ofmaximum Fig. 3. Scatter plot of T= 200 samples from a mixture of Gaussian clusters. Red line: direction of the optimal projection according to the proposed criterion (L1-uLDA). Green line: optimal Fisher’s discriminant direction (LDA). Black line: direction of maximum variance (L2-PCA). Fig. 4. Histogram of the data projected onto the direction determined by L1- uLDA in Fig. 3. The two modes point out the existence of two clusters. variance of the data, constructed by classical L2-PCA, whereas the green line marked ‘LDA’ represents Fisher’s discriminant direction for this dataset. Eventhoughthethreeapproachespresentasimilarbehaviorin this example, strictly speaking classical L2-norm based PCA is not a linear discrimination technique and, therefore, it may yield different results in certain cases. Supporting this claim, Fig. 5 shows the projection directions obtained after repeating the experiment under the same conditions except for the covariance matrices, which are set to V1=V2=diag(1,3), so that the direction of maximum variance lies now along the x2-axis. 3) More Than Two Dimensions: Let us now extend the previous analysis to more than two dimensions. Consider data from two classes with prior probabilities π1=π2=1/2in a p-dimensional space, with p>2. Data from class i,i=1,2, are drawn from a Gaussian distribution N(µi,Vi), with µ1= MARTÍN-CLEMENTE AND ZARZOSO: LDA via L1-PCA OF WHITENED DATA 233 TABLE I ANGLES IN DEGREES BETWEEN LDA AND L1-uLDA DIRECTIONS (ΔθL1), AND LDA AND L2-PCA DIRECTIONS (ΔθL2). VALUES SHOWN REPRESENT MEAN ± STANDARD DEVIATION OVER 10 INDEPENDENT REALIZATIONS OF T= 200pDATA SAMPLES.p:DATA DIMENSIONALITY;μ:DISTANCE BETWEEN CLASS MEANS Fig. 5. Linear discrimination performance for classes with non-isotropic covariance matrices, under the same general conditions of the experiment of Fig. 3. L2-PCA fails to find a discriminant solution, whereas L1-uLDA projector still lies close to LDA’s. −μ 2e1,µ2=μ 2e1, where e1is the vector with 1 in the first entry and 0’s elsewhere. Accordingly, parameter μrepresents the distance between the class means. The covariance matrices are chosen either as (a): V1=V2=I,or(b):V1=V2= diag(1,p,...,p). In both cases, the direction resulting from Fisher’s discriminant is parallel to the x1-axis. In (b), the x1-axis is also the direction of minimum variance of the classes. The LDA projection vector (aLDA), the L2-PCA principal direction of the data (aL2) and the proposed L1-uLDA projection vector (aL1) are estimated from the same T= 200 pdata points and normalized to unit length. To compare Fisher’s discriminant with the two other methods, the angles formed by the LDA vector aLDA with aL1 and aL2 are calculated by the formula Δθx= arccos(|aT LDAax|), where x∈{L1,L2}, for different parameter pairs (p, μ). The absolute value simply forces the angle to lie in the first quadrant. Table I shows the mean angle values ±standard deviations (in degrees) for 10 independent data realizations. Observe that θL1 tends to 0◦as μincreases, which means that aL1 gradually aligns itself with aLDA, achieving a perfect fit when the clusters are well separated, as predicted by our theoretical analysis. While the different approaches indeed become equivalent when the classes are clearly distinct, L1-uLDA keeps closer to LDA than L2-PCA when the classes tend to overlap. This observation is particularly evident in case (b), where aL2 lies almost orthogonal to aLDA for small μ. 4) Data Corrupted With Outliers: Let us return to the case p=10 and μ=5 with V1=V2=Iof the previous experiment. We replace the 20% of the data from each class by outliers generated by raising normalized Gaussian random samples to the third power and centering them on the corresponding class means. Some outliers considerably far from the clusters are obtained in this manner. For the corrupted data, the mean angle between the optimal L1-uLDA projection vector and the horizontal axis, averaged over 10 independent experiments, is equal to 8.8◦with an standard deviation of 4.9◦. This mean angle is equal to 9.4◦for LDA applied on the same corrupted data (standard deviation equal to 3.3◦). Both means are compared with the aid of Student’s t-test and do not differ significantly at the 5% level. For comparison, the angle is 75.9◦(standard deviation 1.4◦) when calculated by L2-PCA. This experiment suggests that the proposed L1-uLDA method is as robust to outliers as the original LDA method while operating in an unsupervised manner. 5) Unbalanced Clusters: The L1-uLDA criterion may fail when one of the classes is much more frequent than the other, so that the cluster masses become too unbalanced. Consider two-dimensional random samples from a mixture of two Gaussian distributions with means µ1=−μ 2e1,µ2=μ 2e1, where e1=[1,0], and common covariance matrix V1=V2=I, where Iis the identity matrix. Unlike the above experiments for classes with the same prior probabilities π1=π2=0.5, now the clusters are not assumed to have a similar number of points. Fig. 6a plots the objective function J(θ), defined in Eqn. (44), for μ=3and different values of the class prior probabilities π1and π2, with π1+π2=1. The plots show that, for unevenly sized clusters (e.g., π1=0.1or π1=0.2), the maxima of J(θ)occur at θ=±90◦, and hence L1-uLDA would produce directions orthogonal to Fisher’s discriminant. This outcome can be interpreted as follows. When the priors πiare highly unbalanced, the zero-mean constraint forces the distribution with larger mass to be concentrated around zero. Therefore, the mean of the projected class is always close to zero as well, no matter on which direction we project the data. This prevents the projected data from forming a cluster where the mean is large relative to the standard deviation, which is required for L1-uLDA to work. Fig. 6b plots the new function J(θ)obtained after setting μ=14, resulting in a clearer cluster separation. Observe that the curve π1=0.2now passes through a maximum at θ=0 ◦, although it is not the global maximum. Therefore,L1-uLDA mayfailagaininthis case.See AppendixC for a further explanation. 6) The Multiclass Case: Consider a p=30–dimensional data set, with c=30classes, and 500 samples per class (implying that the prior probabilities are the same for all classes). The class means µiare generated from a normal distribution with 240 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 68, 2020 [30] R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification. Hoboken, NJ, USA: Wiley, 2001. [31] J. Chen, Z. Zhao, J. Ye, and H. Liu, “Nonlinear adaptive distance metric learning for clustering,” in Proc. 13th Int. Conf. Knowl. Discovery Data Mining, 2007, pp. 123–132. [32] C.H. Q. Ding and T. Li, “Adaptive dimension reduction usingdiscriminant analysis and k-means clustering,” in Proc. 24th Int. Conf. Mach. Learn., 2007, pp. 521–528. [33] S. Wang, J. Lu, X. Gu, H. Du, and J. Yang, “Semi-supervised linear discriminant analysis for dimension reduction and classification,” Pattern Recognit., vol. 57, pp. 179–189, 2016. [34] “MATLAB code for the L1-PCA algorithm proposed in [17],” [Online]. Available: https://sites.google.com/view/miloslab/resources. Accessed on: 2018. [35] R. Fisher, “The use of multiple measures in taxonomic problems,” Ann. Eugenics, vol. 7, no. 2, pp. 179–186, 1936. [36] D. Dheeru and E. Karra Taniskidou, “UCI machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml. [37] “Iris flower data set,” [Online]. Available: https://archive.ics.uci.edu/ml/ datasets/Iris. Accessed on: 2018. [38] “Breast Cancer Wisconsin (diagnostic) data set,” [Online]. Available: https://archive.ics.uci.edu/ml/machine-learning-databases/breastcancer-wisconsin/. Accessed on: 2018. [39] M. Turk and A. Pentland, “Eigenfaces for recognition,” J. Cogn. Neurosci., vol. 3, pp. 71–86, 1991. [40] P. N. Belhumeur, J. P. Hespanha, and D. J. Kriegman, “Eigenfaces vs. fisherfaces: Recognition using class specific linear projection,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 19, no. 7, pp. 711–720, Jul. 1997. [41] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in Proc. IEEE Conf. Comput. Vis. Pattern Recognit., 2016, pp. 770–778 [42] D. E. King, “Dlib-ML: A machine learning toolkit,” J. Mach. Learn. Res., vol. 10, pp. 1755–1758, 2009. [Online]. Available: http://dlib.net [43] G. B. Huang, M. Ramesh, T. Berg, and E. Learned-Miller, “Labeled faces inthewild:Adatabase for studying facerecognitionin unconstrained environments,” [Online]. Available: http://vis-www.cs.umass.edu/lfw/index. html. Accessed on: 2019. [44] N. Dalal and B. Triggs, “Histograms of oriented gradients for human detection,” in Proc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit., Jun. 2005, vol. 1, pp. 886–893. [45] E. Chong and S. Zak, An Introduction to Optimization, 4th ed. Hoboken, NJ, USA: Wiley, 2013. Rubén Martín-Clemente received the M. Eng. degree in telecommunications engineering and the Ph.D. degree with highest distinction in telecommunications engineering from the University of Sevilla, Sevilla, Spain, in 1996 and 2000, respectively. He is currently the Head of the Department of Signal Theory and Communications of the University of Sevilla, Spain. He has been a Visiting Researcher with the University of Regensburg, Regensburg, Germany, in 2001 and 2009 andwiththeUniversityof Nice, France, in 2015, 2016, and 2018, respectively. Among other areas, his research interests include signal processing and machine learningwithemphasis on independent component analysisanditsapplicationto biomedical problems. He has authored or co-authored numerous publications on these topics. He was as a Program Committee Member for several international conferences and was a Program Committee Chair of the 5th International Conference on Independent Component Analysis and Blind Signal Separation in 2004. Vicente Zarzoso received the graduate degree with highest distinction in telecommunications engineering from the Polytechnic University of Valencia, Valencia, Spain, in 1996. He began the Ph.D. studies with the University of Strathclyde, Glasgow, U.K. In 1999, he received the Ph.D. degree from the University of Liverpool, Liverpool, U.K. He received the Habilitation to Lead Researches (HDR) from the University of Nice Sophia Antipolis (now member of the Université Côte d’Azur, UCA), France, in 2009. From 2000 to 2005, he held a research fellowship awarded by the Royal Academy of Engineering, U.K. Since 2005, he has been withthe Computer Science, Signalsand Systems Laboratory of SophiaAntipolis (I3S), UCA, CNRS, France, where he is a Full Professor and since 2016 he has been the Head of the “Signals, Images and Systems” (SIS) research team. His research interests lie in the areas of signal processing and machine learning with emphasis on matrix and tensor factorizations, principal/independent component analysis and related techniques, including theoretical aspects and applications in biomedical problems and digital communications. He was an Associate Editor for the IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS from 2011 to 2015 and has been a Program Committee Member of several international conferences. He was a Program Committee Chair of the 9th International Conference on Latent Variable Analysis and Signal Separation (LVA/ICA-2010) and a keynote lecturer at the LVA/ICA-2015 Summer School. He is a Senior Member of the IEEE and an Honorary Member of the Institut Universitaire de France.