Dimensionality reduction plays an important role in pattern classification and has been successfully applied to many realworld applications, such as social media, business, computer vision, bioinformatics, etc. [1]. Among the most widely used supervised dimensionality reduction methods is the Linear Discriminant Analysis (LDA) that aims to find the optimal projection by simultaneously maximizing the between-class variance and minimizing the within-class variance in the projected subspace. The algorithm of solving the traditional LDA objective is to optimize a trace ratio problem that can be transformed into a tractable ratio trace form [2]. This transformation allows us to easily find a closed-form solution by using the generalized eigenvalue decomposition method [3]. Unfortunately, this solution often deviates from the original objectives when outlier samples are present, which can result in uncertainties within subsequent classification procedures [2]. Another problem with this method is that inverting a singular within-class variance matrix can cause the LDA formulation to be ill-conditioned, and in fact this is common when the training data are sampled from highdimensional space or there exist strong-correlated features. The efforts to overcome these issues by optimizing the ratio trace transformation cannot thoroughly solve this singularity problem. Motivated to tackle the original trace ratio problem, many existing works [4]- [6] studied how to directly optimize the trace-ratio objective function.
One major limitation of the traditional LDA is its dependency upon the squared 2 -norm ( 22 ) distance that is highly vulnerable to the presence of outliers or noisy data. This is because the least square criterion can remarkably enlarge the effect of noise on the total loss. Unfortunately, noise appears more or less everywhere in real-world applications. Therefore, it is critical to improve the robustness of the model when random noise or even adversarial attack are present in the data. In contrast to regular noises that are usually small in magnitude and sometimes even invisible to human eyes, adversarial attacks are intentionally constructed by injecting small perturbations into the original data such that the model can be misled to output discouraged results [7]. A more recently used adversarial attacks for modeling real-world data are the adversarial patch attack that couples the original data with a human-manipulated patch [8]. The effectiveness of such attacks can be explained in terms of the black-box principle where the attacker knows nothing about the configurations of the studied model.
To enhance the robustness against noise, many LDA objectives changed to use other distance functions rather than the 2 2 -norm distance. Wang et al. [9] proposed to measure the two variance matrices using a rotational invariant 1norm. However, the runtime of this method is as high as the intractable memory cost of storing the representations because the greedy learning strategy proceeds one projection vector at a time. Tao et al. [10] proposed to select the discriminative patterns by using a 2,1 -norm regularization that is a row-sparsity constraint on the transformation matrix. Since this method fails to uncover the most important features in the given tasks and the learned representations are largely governed by the selected number of dimensions. Wen et al. [11] introduced an orthogonal matrix to connect the original features and transformed features, in which the main information of the original data can persevere in the discriminant subspace. However, their optimization methods easily fall into locally optimal solutions. Zheng et al. [12] tried to build a 1norm discriminant analysis (L1-LDA) using Bayes error bound optimization, but the implementation of the Bayes error bound is not easy as it requires each class to contain the same number of samples.
Recently, there are also several efforts that focus on how to develop robust lost functions [13]- [15]. Although these methods are able to improve the robustness against noisy data, their implementations are usually limited to two notorious problems. One is that some optimization algorithms just simply guarantee the decrease of the objective value, but not the algorithm convergence by nature. Another problem is that the derived solutions using their methods can easily get stuck in local minima rather than landing exactly on the global minima. Moreover, these loss functions were developed with their particular properties, each of which has its own independent variables to adjust the robustness. Hence, in this paper, we aim to provide a novel and robust objective for LDA that has better generalization with a shared tractable and efficient optimization algorithm.
We also noticed that most existing robust LDA objectives do not tackle the mean calculation in a correct way. In many cases, the mean of the data is not zero, and using 2 2 -based distance to calculate the mean value is incorrect when the objective depends on other norm based distance [16]. In this paper, we thus do not only propose a novel generative LDA model that integrates a very broad family of robust loss functions, but also embed the mean calculation into the objective function that can remove the optimal mean automatically during learning time. Since existing optimization algorithms are not able to well handle our generative objective with the non-smooth and non-convex property, as one important algorithmic contribution, we further derive an efficient solution algorithm that can guarantee both the objective and sequence convergences. Through our mathematical analysis, our proposed algorithm is proved to find a globally optimal solution at a sub-linear rate. In our extensive experiments, our proposed method outperforms several competitive dimensionality reduction approaches and has a great advantage of handling various noisy date or adversarial attacks.
Given the training dataset of X ∈ R d×n , we can represent it using a few sub-design matrices X = [X 1 , ..., X c ],
where
) corresponds to a collection of data points belonging to the ith class. The traditional binary-class LDA method assumes that, given data points X i = [x i 1 , ..., x i ni ] in the i-th class and data points X j = [x j 1 , ..., x j nj ] in the j-th class, both multivariate probability densities are multivariate Gaussian distributions with arbitrary mean vectors μ and a covariance matrix Σ, that is:
, and
, where μ i is the mean value of data points in class i, and μ is the mean value of the whole data points. Due to the homogeneity assumption that Σ i = Σ j = Σ xx , the ratio of the two densities can be written as:
.
(
Since the logarithms is a monotonically increasing function, Eq. ( 1) follows that:
where μ = μi+μj 2
. Also, assume that given a randomly selected data point x belonging to the i-th class, we have the prior probability:P (x ∈ X i ) = π i . Using the Baye's theorem, we have the logarithms of the ratio of the posterior probabilities:
where
). Since the function L(x) is a linear function of x, the classes i and j can be separated from the linear combination, w T X, of the data matrix X. In other words, the objective of binary-class Linear discriminant analysis (LDA) is to find a projection w that can maximize the squared distance of the two class means, w ∝ Σ -1 xx (μ iμ j ), comparing to the within-class variation of that distance. From this point of view, in multiclass LDA, we define the between-class scatter matrix as
, where n i is the total number of data points in class i, μ i is the mean value of data points in class i, and μ is the mean value of the whole data points. The objective of multiclass Linear Discriminant Analysis is to find a transformation matrix W ∈ R d×m that can minimize the squared distance of data points within the same class and maximize the squared distance of data points between different classes. We can thus represent the objective function of multiclass LDA as:
(4) Rather than measuring the between-class scatter matrix S w , we also need to measure the total scatter matrix S t = n i=1 (x iμ)(x iμ) T , given the fact that S t = S w + S b . The objective function in Eq. ( 4) can be rewritten as:
However, this objective is challenging to optimize due to the fact that the numerator is a minimization problem while the denominator is a problem to maximize variance in the projected space, max
A smart way of solving this challenge is to convert the maximization problem into a minimization problem by using the reconstruction errors [6], [17]. We can therefore change it into:
The 2 2 -norm used in this objective is known as the Euclidean norm. In many contexts, the 2 2 is undesirable because of the great sensitivity to noisy data or outliers [18]. In addition, as we discussed in the linear model L(x) in Eq. ( 3), this squared error is motivated by the data with a Gaussian prior. However, [19] proposed that if the linear model changes to build upon a Laplace prior, this method can provide a more robust algorithm and produce lower error rates than previous techniques with a Gaussian prior. With this recognition, it is easy to see that the problem in Eq. ( 6) with a Laplace prior distribution reduces to the 1 -norm optimization problem. We also note that the fragile 2 2 -norm distance can be replaced with multiple robust loss functions (for example, Nie adaptive loss function [13], Log-cosh loss function, Huber loss function [20], Cauchy loss function or Lorentzian function [21], or Barron adaptive loss function [14], etc., as shown in Table I, we propose to combine the LDA objective with such robust loss functions. In the most practical scenarios, optimization algorithms used for solving such robust loss functions confront several intractable challenges as discussed in the introduction section. We thus aim to present a generalized LDA formulation that consists of different loss functions, and meanwhile, this objective can be solved using shared and tractable optimization algorithm. From this point of view, we define a generative loss function L(•) that can be replaced by a family of loss functions like 2 , 1 , Huber, etc, as listed in out Table 1, all of which is required to fulfill an important property -Lipschitz continuous gradient. Following intuitions above, Eq. ( 6) can be changed into:
Likewise, the numerator in Eq. ( 5) can be changed into:
We can now view Eq. ( 7) as a constrained condition to Eq. (8). Specifically, if L = I, this constraint changes back to the Euclidean norm. If L = 0, all data points in this constraint collapse into a single point, which is useless in realworld scenarios. To avoid this collapse problem, we change the formulation to an inequality constraint:
where c is a positive constant. To obtain a more convenient but equivalent problem, we observe that taking any positive values of c just results in the L being replaced by c 2 L, meaning that the inequality constraint in Eq. ( 9) can be simplified into an equality constraint [22]. We can now write our optimization problem as:
This constrained minimization problem can be solved by constructing a Lagrange function, which is given by:
This leads to the optimization problem as below:
where the new variable λ is a multiplier for the constraint.
It is worth noting that the terms μ and μ i are often obtained by μ i = ni j=1 x i j and μ = n i=1 x i . However, the intuition behind this calculation depends on the Euclidean distance, and, as a result, can be invalid when the loss function is built upon some robust norm based distance, such as 1 -norm or 2,1norm, according to [16]. This is because that the mean of a set of input data depends on the definition of distance. Different distance metrics/loss functions lead to different mean. Also, using different loss functions is equivalent to use to a different distance. At this point, the mean in our method is a variable, but not a fixed value like the traditional one which is defined for the 2 -norm distance, since our objective is designed for representing a family of different robust loss functions. We thus propose a novel generative LDA objective function that can automatically remove the optimal mean from the data over the iterations.
Learning Eq. ( 12) can be difficult in practice as the objective is non-convex and non-smooth, which cannot admit closed form solutions. Rather than minimizing L with respect to the entire parameters simultaneously at each time step, we can think of minimizing L alternatively with respect to a single parameter at a time. For example, if we minimize L with respect to the parameter W , then minimize it with respect to μ and so on, repeatedly cycling through all parameters, we are guaranteed to obtain an optimal minimum. This practice is called as Alternating Minimization method (AM). As a result, when applying the AM to the objective in Eq. ( 12), we obtain the solution update rule as:
where k is the time step index, and W, M i , M, Λ are the sequences of associated variables.
Formula Gradient
Huber loss
L(y, ŷ) = log ( y-ŷ c 0
Log-cosh loss L(y, ŷ) = log cosh(y -ŷ)) ∇L = tanh(y -ŷ) = e 2(y-ŷ) -1 e 2(y-ŷ) +1
Nie adaptive loss L(y, ŷ) =
otherwise.
The matrix factorization in Eq. ( 13) is still a non-convex problem due to the constraint on each sequence set. An efficient way that such a constraint can be used in the service of regularizing the optimization function is by adding a proximal regularization term to the update. This method has been used primarily in the context of nonconvex and nonsmooth problems [23]- [25]. Using the proximal terms introduced in the alternating linearized minimisation, we denote our regularized optimization problems as:
where ∇L(W k ) means the gradient of Eq. ( 11) at the point of W (k) , and < ∇L(W (k) ), W -W (k) > means the dot product between the gradient ∇L(W (k) ) and the offset (W -W (k) ).
Another important advantage of this method of adding the proximal term to optimization problem is that the regularization terms can avoid the updated solution being changed to a value that differs significantly from the value of previous solutions.
We now discuss the solution algorithm in more details, with emphasis on its application to a variety of loss functions defined in the generative loss function. Table I shows a set of loss functions that can be incorporated into the generative loss function, including, but not limited to the 2 2 -norm, Huber loss function [20], Nie adaptive loss function [13], Barron adaptive loss function [14], etc. Many extensions to this generative loss function are possible, as long as they remain the Lipschitz continuous gradient. To give a sense of how to update the solution given in Eq. ( 14), we need to define two notations by simplifying
The gradient of A i and B with respect to a variable is written as grad(A i ) and grad(B), each of which is derived using the element-wise product. For example, suppose we want to train our objective that incorporates Barron adaptive loss function [14] with the use of λ = 0, the general gradient listed in Table I (defined as grad(Y ) in this case) can be written as grad
, where is the Hadamard division, • is the Hadamard product.
Step 1: Updating W (k+1)
where
We can easily obtain this solution by finding the nearest orthogonal matrix to the given matrix M (k) , known as the orthogonal Procrustes problem [26],
where P 1 and P 2 are from the singular value decomposition
where
. Setting its derivative equal to 0 we have the solution:
Step 3: Updating μ (k+1)
where
Setting its derivative equal to 0 we have the solution:
where
The solution can be obtained by setting its derivative equal to 0:
Algorithm 1: Proximal Alternating Linearized Minimization The solution algorithm to solve the generative LDA objective in Eq. ( 10) by using the proximal alternating linearized minimisation is summarized in Algorithm 1. It is worth noting that our new solution algorithm guarantees that both the objective and the solution sequences converge to the globally optimal solution at a sub-linear convergence rate. Due to space limitation, the proof of the convergence of Algorithm 1 will be provided in the extended journal version of this paper.
We evaluated the performance our model on several real data benchmarks including COIL20, FERET, USPS, MNIST, and Olivetti. [27]- [29]. As our model aims to improve the model robustness against noise, we were motivated to design a few experiments with noisy data, in which block disturbance or Gaussian noise was added into the input data. Several state-of-the-art methods were also implemented for comparison, involving RSLDA [11], 2,1 -LDA [30], GLDA [31], LDA, stacked Restricted Boltzmann machine (RBM) [32], and Stacked AutoEncoder (SAE) network [33].
Datasets descriptions. The COIL20 dataset contains 1440 images of 20 subjects, in which each subject has 72 images taken at pose intervals of 5 angle degree. The size of each image is 128 × 128. The FERET dataset contains more than 14, 000 facial images of 1199 individuals. Since the main purpose of our experiment is to test the object recognition against noise attacks, we randomly selected a subset of the FERET dataset for our experiment. This gives rise to 1400 images of 200 individuals, each of which has 7 face images. The size of the image was converted to 60 × 60 pixels. The USPS image dataset is commonly used for handwritten text recognition research. This dataset contains a total of 9, 298 samples and each sample has 16 × 16 greyscale pixels. The MNIST dataset is another large database of handwritten digits that contains 70, 000 samples. Each greyscale image is 28×28, representing the digits 0 -9. The Olivetti (or AT&T) dataset consists of 400 different images of 40 distinct individuals. Each individual has 10 different images taken at different conditions, and each image was resized to 64 × 64 pixels.
Experimental results and analysis. We compared our model against seven other methods involving RSLDA [11],
2,1 -LDA [30], GLDA [31], LDA, stacked Restricted Boltzmann machine (RBM) [32], and Stacked AutoEncoder (SAE) network [33]. In this experiment, we demonstrate the robustness of our model on different noisy data sets. The loss function we used in our model is Nie and Barron adaptive loss function, respectively. Our training procedures were implemented using the stratified 6-fold cross-validation that partitions data into k non-overlapping folds to fit the train/test set by preserving the percentages of samples for each class. We repeated the stratified 6-fold cross-validation six times and then reported the mean performance across all folds and all repeats. See Table II for a comparison of the performance of our method and seven other models on noisy datasets. Note that the noisy data used in the our tables were formed by adding Gaussian noise.In these cases, the Gaussian noise was randomly distributed with the variance of 0.05. The dimensionality of the above data was reduced to a fixed level that is 19, 25, 9, 9, and 25 for COIL20, FERET, USPS, MNIST, and Olivetti datasets, respectively. The table report the average classification accuracy, precision, and recall scores over 50 runs with their standard deviations. We can see that our model clearly outperforms other competitive methods.
We propose a generative LDA learning method in which multiple robust loss functions can be interchangeable. Also, to solve the singularity problem in a natural way, our objective depends on the trace-ratio formulation that can avoid the inversion of the within-class variance matrix. In addition, our objective automatically center the data by quantifying the optimal mean value during the learning process. We further introduce an effective algorithm to solve the resulting nonconvex and non-smooth problem. Through our mathematical analysis, we conclude that our optimization algorithm can find globally optimal solutions with guaranteed convergences of both the objective value and the solution sequence, and the converging rate is as fast as sub-linearity. The experimental results on several real-world datasets shows that our model consistently outperformed its state of the art counterparts.