Robust Kernel Dictionary Learning Using a Whole Sequence Convergent Algorithm

1 Introduction

Sparse coding, which assumes that the signal can be sparsely represented by a dictionary, has become an active topic for scholars in the fields of machine learning and signal processing[Cheng et al., 2009]. However, the linear reconstruction assumption which is explicitly imposed by many existing work is not valid for many practical signals[Cheng et al., 2013][Yang et al., 2015], particularly those well-modeled by manifolds[Hussein et al., 2013][Ma et al., 2014]. [Gao et al., 2010] showed that one can use the kernel trick to perform sparse coding in the high-dimensional feature space instead of the input space in order to capture the non-linear structure of signals more efficiently. Kernel sparse coding has therefore extensive applications in many fields. On one hand, the high dimensionality of the feature space typically yields a more discriminative representation than input data space[Gao et al.,

∗This work was supported in part by the National Key Project for Basic Research of China under Grant 2013CB329403; in part by the National Natural Science Foundation of China under Grants 61210013 and 61273256; and in part by the Tsinghua University Initiative Scientific Research Program under Grant 20131089295.

2013][Nguyen et al., 2013]. On the other hand, kernel trick provides effective strategy for sparse coding on Riemannian manifold[Harandi and Salzmann, 2015], of which structure is more complicated than the Euclidean space[Wang et al., 2014].

Similar to the linear sparse coding, the dictionary learning problem of kernel sparse coding is also very important. If we map the training samples into a higher dimensional feature space by a proper mapping function, and learn the dictionary in this feature space, we expect a better linear representation. In [Gao et al., 2010], the dictionary learning problem was tackled by gradient descent over the basis set in the input space directly. To compute the gradient, they had to resort to a specific kernel function such as the Gaussian kernel. Similarly, [Harandi et al., 2012] and [Li et al., 2013] addressed the dictionary learning problem on Symmetry-Positive-Definite (SPD) manifold, using Stein kernel and Log-Euclidean kernel, respectively. Very recently, [Zhang et al., 2015] investigated the online dictionary learning on SPD manifold. All of those work deal with specified kernel and can limit their applicability considerably. In [Xie et al., 2013], a nonlinear generalization of sparse coding and dictionray learning on manifold was proposed. This method is not based kernel and requires an extra affine constraint on the coding coefficients, which may be unexpected in practice. [Li and Ngom, 2012] revealed that the dictionary optimization only needs inner products of samples and this property can be used to kernelize linear sparse coding. Recently, some scholars proposed more principled kernel dictionary learning method[Anaraki and Hughes, 2013][Nguyen et al., 2013][Liu et al., 2014][Kim, 2014][Harandi and Salzmann, 2015]. Those methods are based on the basic conclusion that the dictionary atoms in the feature space can be linearly reconstructed by the samples in the feature space. As a result, the complicated dictionary learning problem in feature space is formulated as the rather simpler search of a coefficient matrix.

Nevertheless, existing kernel sparse coding and dictionary learning methods utilize square loss to measure the reconstruction capability. Such a loss function may be sensitive to the outliers[Liu et al., 2015]. Although the above literatures proposed various methods for dictionary learning, none investigated the outliers which may lie in the training sample set. It should be noted that in [Nguyen et al., 2013], the authors concerned the robustness but only added the corruptions on the testing samples and the training samples are still clean. Therefore, the robust kernel dictionary learning problem, which aims to learn a dictionary in the feature space while isolating the outliers, has not been addressed. As a comparison, the robust dictionary learning problem on Euclidean space has been extensively studied by many scholars. In [Chen and Wu, 2013], the error source decomposition technology was developed to solve the robust dictionary learning, and [Pan et al., 2014] investigated the robust non-negative dictionaries learning. It is not clear how to extend their work to the kernel case. In [Nie et al., 2013], a robust kernel dictionary selection method was developed for active learning. However, this work addressed the problem of dictionary selection, but not dictionary learning. In addition, [Kong et al., 2011] addressed the robust nonnegative matrix factorization. [Xia et al., 2012] investigated the robust kernel nonnegative matrix factorization. Both work did not consider the influence of sparsity.

In this work, we address the robust dictionary learning problem under the framework of kernel sparse coding. The main contributions are listed as follows:

1. A new optimization model is proposed to learn the robust kernel dictionary while isolating outliers in the training samples. This model is essentially based on the decomposition of the reconstruction error into small dense noises and large sparse outliers. The latter is formulated as the product of the sample matrix in the feature space and a diagonal coefficient matrix. This facilitates learning kernel dictionary.
2. A whole sequence convergent algorithm is developed to solve the non-convex robust kernel dictionary learning problem. The solution sequence is proved to be a Cauchy sequence and convergences to the critical point of the original optimization problem.
3. We perform empirical comparisons of the proposed method with the existing dictionary learning methods in the applications to image classification, which justifies that our method yields more robust results.

The rest of the organization of this paper is as follows: Section 2 gives the problem formulation, followed by the algorithm design and analysis in Section 3. Section 4 gives the experimental results.

Notations: For a matrix M. The matrix norm ||M||2 is defined as ||M||2 = √T r(MT M). The pseudo-norm ||M||0 is defined as the number of the non-zero element in M. The pseudo-norm ||M||0,2 is defined as the number of the nonzero columns in M. For a vector u, we use ||u||2 to denote its 2-norm, and ||u||0 to count the number of the non-zero elements. The symbol I represents the identity matrix with compatible dimension.

2 Problem formulation

Given the N training samples {yi}Ni=1 ⊂ M, where M is a Riemannian manifold which also can be a subset of a Euclidean space. If we map the training samples into a higher dimensional space by a proper mapping function, we expect a better linear representation. The linearity in feature space corresponds to the nonlinearity in the original space. To this end, we denote Φ(·) : M→ H to be the implicit nonlinear mapping fromM into a high-dimensional (maybe infinite dimensional) dot product space H. This mapping function is associated with some kernel κ(yi, yj) = ΦT (yi)Φ(yj), where yi, yj ∈ M. For convenience, we denote the dimension of H as d, which may be infinite.

The aim of dictionary learning is to empirically learn a dictionary adapted to the training sample set {yi}Ni=1 that we want to sparsely code. Therefore we need to determine some atoms d1, d2, · · · , dK ∈ H, where K < N is the size of the dictionary, to sparsely represent each training sample in the feature space. By denoting Φ(Y) = [Φ(y1), ...,Φ(yN)] ∈ Rd×K and D = [d1, ..., dK] ∈ Rd×K , we can formulate the kernel dictionary learning problem as

where C ∈ RK×N is the sparse coding matrix and λ1 is the sparsity penalty parameter.

By using the mapping function Φ(·), we can transform the problem on Riemannian manifold to the linear coding problem in feature space. This is the great advantage of the kernel sparse coding[Gao et al., 2013][Harandi and Salzmann, 2015]. While achieving great success in nonlinear sparse coding, such a formulation admits challenge to the dictionary learning since the dictionary atoms dj may be in infinite dimensional space. Fortunately, some recent literature pointed that the dictionary can be represented by D = Φ(Y)A, where A ∈ RN×K is a coefficient matrix. This means that the dictionary atoms can be linearly reconstructed by the training samples in the feature space. This conclusion was proved in [Nguyen et al., 2013] and [Kim, 2014]. Based on this formulation, the dictionary learning problem becomes

To avoid the non-unique solution, we impose the constraint that ||aj||2 = 1, where aj is the j-th column of A. Such a formation provides significant convenience since the learning of dictionary becomes the search of the matrix A and provides a principled derivation for nonlinear dictionary learning and sparse coding that essentially reduces to linear problems for any type of kernel function.

However, the reconstruction error term in Eq.(2) is the squared loss function in the feature space and therefore it is very sensitive to the data outliers, which may dominate the objective function. To attenuate the influence of outlier, we introduce an error matrix E = [e1, · · · , eN ] ∈ Rd×N , where ei ∈ H for i = 1, 2, · · · ,N, to give

where λ2 is the parameter to penalize the column sparsity of E. The motivation of this model is to isolate some few columns in Φ(Y), which cannot be well reconstructed by the dictionary. However, the challenge in solving model (3) lies in the fact that ei is in the feature space and may be infinite dimension. To this end, we try to represent E using the samples Φ(Y). Since the non-zero columns of E correspond to the columns of Φ(Y), we can represent

where R ∈ RN×N is a diagonal matrix. As a result, the problem of outlier search is transformed into the problem of determination of R. Furthermore, we assume that there does not exist all-zero column in Φ(Y), i.e., ||Φ(Y)||0,2 = N, then we have

where r = diag(R) is defined as a column vector accommodating the diagonal elements of the diagonal matrix R.

Therefore the robust kernel dictionary learning problem can be formulated as the following optimization problem:

Please note that although the number of outlier is small, we cannot impose the sparsity on R since it is a diagonal matrix and therefore it is highly sparse. Instead, we encourage the sparsity of the vector which is formed by the diagonal elements. This shows that the outlier isolation term in (4), key to the robust dictionary learning, can be kernelized. Please note that the reconstruction error term now can be formulated as ||Φ(Y)(I-AC-R)||22 = Tr{(I-AC-R)T KY Y (I-

The model in (4) provides twofold roles: (1) Learn the dictionary coefficient A from the training sample set; and (2) Isolate the outliers. However, the model (4) is highly nonconvex due to the coupling term AC and the zero norm terms imposed on C and r. Neglecting the robustness problem, [Harandi and Salzmann, 2015] and [Nguyen et al., 2013] developed method for dictionary learning. Both method are essentially based on the Method of Optimal Directions (MOD) method and Kernelized KSVD (K-KSVD) method. However, no whole sequence convergency can be guaranteed for such algorithms. That is to say, those algorithms at most can be proved that the functional value decreases at each iteration, but the solution sequence itself may not convergent. Motivated by the alternating proximal linearized method proposed by [Bolt et al., 2013] and [Bao et al., 2014], we develop a whole sequence convergent algorithm to solve the optimization problem in (4).

3 Algorithm design and analysis

To design the whole sequence convergent algorithm, we de-

where IA(A) denotes the indicator function of A which satisfies IA(A) = 0 if A ∈ A and +∞ otherwise, where A = {A ∈ RN×K : ||aj||2 = 1, j = 1, 2, · · · ,K}. Similarly, the symbol IC(C) denotes the indicator function of C, where C = {C ∈ RK×N : |Cij| ≤ cm}. The symbol IR(R) denotes the indicator function, where R = {R ∈ RN×N : |Rii| ≤ rm for i, j = 1, 2, · · · ,N, and Rij = 0 for i ̸= j}. The values of cm and rm are simply set to serve as the upper bounds.

Motivated by the alternating proximal linearized method[Bao et al., 2014], we need to alternatively update C, A and R. Before further processing, we formally define the proximal operator as

In the following we introduce how to solve the three unknown matrices C, A and R. For conveniens, we use the superscript (k) to denote the values at the k-th iteration.

3.1 Calculating C

This step essentially realizes sparse coding. At the (k + 1) iteration, we are given C(k) , A(k) , R(k) , then we have

where µ(k) c is appropriately selected step size, and the gradient is calculated as

The (i, j)-element of C(k+1) can then be easily obtained as

3.2 Calculating A

This step tries to update the dictionary coefficient matrix. Since A is column normalization, we can update it column by column as a(k+1) j ∈ ProxF( ¯A(k) ) µ(k) aj (a(k) j -

In addition, µ(k) aj is appropriately selected step-size and the gradient is calculated as ∇aj H(A,C,R) = -2KY Y (I - R)CT qj + 2KY YACCT qj, where qj is a K-dimensional vector of which the j-th element is 1 and the other elements are zeros.

we have

and the solution is easily obtained as a(k+1) j = a∗j /||a∗j ||2. After obtaining a(k+1) j for j = 1, 2, · · · ,K, we can naturally obtain the dictionary coefficient matrix solution A(k+1) = [a(k+1) 1 , · · · , a(k+1) K ].

3.3 Calculating R

Given C(k+1) , A(k+1) and R(k) , the value of R can be updated as R(k+1) ∈ ProxQµ(k) r (R(k) - 1 µ(k) r ∇RH(A(k+1) ,C(k+1) ,R(k) )), where µ(k) r is appropriately selected step size and the gradient is calculated as

Denote R∗ = R(k) - 1 µ(k) r ∇RH(A(k+1) , C(k+1) , R(k) ), then we have

The solution can be obtained as

3.4 Selection of step sizes

The step-sizes µ(k) c , µ(k) aj , µ(k) r should be determined during the solving procedure. To this end, we need to calculate the Lipschitz constants L(k) c , L(k) aj and L(k) r which satisfy ||∇CH(A(k) ,C1,R(k) ) - ∇CH(A(k) ,C2,R(k) )||2 ≤ L(k) c ||C1 - C2||2, ||∇aj H(aj1,C(k) ,R(k) ) - ∇aj H(aj2,C(k) ,R(k) )||2 ≤ L(k) aj ||a(k) j1 - a(k) j2 ||2, and ||∇RH(A(k) ,C(k) ,R1) - ∇RH(A(k) ,C(k) ,R2)||2 ≤ L(k) r ||R1 - R2||2, respectively. By simple calculation, we can select L(k) c = ||A(k)T KY YA(k) ||2, L(k) aj = 2||KY Y ||2||C(k)C(k)T qj||2, and L(k) r = 2||KY Y ||2. After obtaining the values of L(k) c , L(k) aj , and L(k) r , we can determine the step-sizes as µ(k) c =max(ρL(k) c , µmin), µ(k) aj = max(ρL(k) aj , µmin) and µ(k) r =max(ρL(k) r , µmin), where ρ > 1 and µmin is a prescribed upper bound. In our work, we set

ρ = 1.1 and µmin = 0.1.

From the above derivation we can clearly see that the sequences {L(k) c }, {L(k) aj } and {L(k) r } are all bounded sequence because the sequences {C(k)} and {A(k) }, and the value of KY Y are bounded.

3.5 Algorithm analysis

The iteration procedure developed in the above sections presents excellent property which is summarized in the following theorem.

Theorem 1: The solution sequence {A(k) ,C(k) ,R(k) }, which is generated by the iteration procedure in Eqs.(7), (8) and (10), is a Cauchy sequence and converges to a critical point of (4).

Proof: The proof is essentially based on the results in [Bao et al., 2014] and [Bolt et al., 2013]. Here we provide a straightforward sketch proof. First, the objective functions H(A,C,R), F(A), G(C) and Q(R) are obviously semialgebraic functions and therefore the whole objective function in (4), which can be written as H(A, C, R) + F (A) + G(C) + Q(R), is a semi-algebraic function, and therefore satisfy the Kurdyka-Lojasiewicz property. Secondly, the sequence {A(k) ,C(k) ,R(k)} is bounded and the step-sizes µ(k) c , µ(k) aj and µ(k) r are bounded. Thirdly, ∇H(A,C,R) = (-2KY Y (I-R)CT +2KY Y ACCT , -2AT KY Y (I-R)+ AT KY Y AC, -2KY Y + 2KY Y AC + 2KY Y R) has Lipschitz constant on any bounded set. Therefore, by adopting the results of Theorem 6.1 in [Bao et al., 2014], the solution sequence is a Cauchy sequence and converges to a critical point of (4).

The conclusion in Theorem 1 clearly illustrates that the sequence {A(k) ,C(k) ,R(k) } is a sequence whose elements become arbitrarily close to each other as the sequence progresses. Therefore, the developed algorithm achieves whole sequence convergence. Such a property provides important practical values since the number of iterations does not need to be determined empirically[Bao et al., 2014].

4 Experimental results

In this section we present several experimental results demonstrating the effectiveness of the proposed robust kernel dictionary learning.

4.1 Synthetic Data

The synthetic data is used to illustrate the outlier rejection role of the proposed method. Similar to the setting in [Nguyen et al., 2013], we learn the dictionary from 500 data sample set {yi}500 i=1 that is generated from a 2-dimensional parabola {yi = [yi,1, yi,2]T ∈ R2 |yi,2 = y2 i,1}. To corrupt the data, we replace yj for j = {120, 220, 320, 420} as newly generated random vectors. Therefore 4 outliers are incorporated into the sample set. We then preprocess all the samples to have unit L2 norm,

We use a polynomial kernel of degree 2 and set the size of the dictionary to be 20. For K-KSVD, the sparsity is set to 5. For the proposed method, we solve the optimization problem in Eq.(2) with λ1 = 0.01 and λ2 = 0.001.

Fig.1 shows the results. The left panel shows the convergence behavior of the proposed method. It shows that the increment indeed converges to zero after several iterations. We also show the convergence behavior of the sparse coding coefficient matrix of K-KSVD method in the right-top panel, which illustrates that the generated sequence does not converge to zero. This validates that K-KSVD has at most subsequence convergence, while the proposed method has whole sequence convergence. Finally, we show the elements of the vector r which is the diagonal elements of R in the rightbottom panel. It clearly shows that the imposed outlier can be reliably isolated, while K-KSVD has no capability to isolate outliers. The results clearly indicate the effectiveness of the proposed algorithm.

4.2 Dataset

In this section we use 4 publicly available datasets to compare several dictionary learning methods. For USPS dataset, we use the vector feature and Gaussian kernel κ(yi, yj) = exp(-γ||yi - yj||22 ). For the rest three datasets, we calculate specified feature vector fu,v for each pixel at (u, v) in one image and use the popular Region Covariance Descriptor (RCovD)[Hussein et al., 2013] to represent the image. Since RCovDs do not lie in Euclidean space, we use the Log-Euclidean kernel κ(yi, yj) = exp(-γ||logm(yi) - logm(yj)||22 )[Li et al., 2013]. By such settings we show that the kernel dictionary learning can be adapted different kernel forms. In all of the experiments, the parameter is set as γ = 0.02. In the following we list some details about the

datasets.

1. USPS digit recognition[Nguyen et al., 2013]: This dataset contains 10 classes of 256-dimensional handwritten digits. For each class, we randomly select N = 500 samples for training and 200 samples for testing. This setting was suggested in [Nguyen et al., 2013] for K-KSVD method.
2. Virus image classification[Harandi and Salzmann, 2015]: This dataset includes 15 different classes. Each class has 100 images of size 41 × 41. For each class, we randomly select N = 80 samples for training and 20 samples for testing. At each pixel (u, v) of an image, we compute the 25-dimensional feature vector fu,v = [Iu,v, |∂I/∂u|, |∂I/∂v|, |∂ 2 I/∂u2 |, |∂ 2 I/∂v2 |, |G0,0 u,v |, · · · , |G4,5 u,v |]T , where Iu,v is the intensity value, Go,s u,v is the response of a 2D Gabor wavelet with orientation o and scale s, and | · | denotes the magnitude of a complex value. Therefore, we generated 20 Gabor filters at 4 orientations and 5 scales.
3. Kylberg texture classification[Kylberg, 2011]: This dataset includes 28 classes. Each class has 160 samples which are resized as 128 × 128. For each class, we randomly select N = 80 samples for training and 80 samples for testing. At each pixel (u, v) of an image,
4. we compute the 5-dimensional feature vector fu,v = [Iu,v, |∂I/∂u|, |∂I/∂v|, |∂ 2 I/∂u2 |, |∂ 2 I/∂v2 |]T , where Iu,v is the intensity value.
5. UCMERCED scene classification[Yang and Newsam, 2010]: This dataset includes 21 challenging scene categories with 100 samples per class. For each class, we randomly select N = 80 samples for training and 20 samples for testing. At each pixel (u, v) of an image, we compute the 15-dimensional feature vector fu,v = [f TR,u,v , f TG,u,v , f TB,u,v ]T , where f T C,u,v = [IC,u,v, |∂IC/∂u|, |∂IC/∂v|, |∂ 2 IC/∂u2 |, |∂ 2 IC/∂v2 |] and IC is the intensity image for the C channel and C ∈ {R, G, B} represents one of the color channel.

In all of the above datasets, the training/testing split is randomly repeated for 10 times and the average results are reported. In addition, we generate corresponding corrupt datasets by selecting ρ% training samples in each class and replace them with the corresponding outlier images. The outlier image is produced by overlaying a rand noise block of which area is a half of the original image. The value of ρ is selected from {0, 10, 20, 30, 40, 50}.

For comparison, we designed the following dictionary learning methods:

1. randDict: This method just randomly selects K atoms subset {s1, · · · , sK} from the training sample set {yi}Ni=1 to construct the dictionary.
2. K-Means method: For vector descriptor, we just use the conventional k-means clustering method to get the K dictionary atoms. For RCovD, we first map the covariance matrices to the linear space by matrix logarithm, in which the clustering is performed and the results are then mapped back to the Symmetry-PositiveDefinite manifold. As a result, we also get the dictionary as {s1, · · · , sK}.
3. GRADient method. For the dataset using Gaussian kernel, we follow the method in [Gao et al., 2010] to learn the dictionary. For the datasets using Log-Euclidean kernel, we follow the method in [Li et al., 2013] to learn the dictionary. The dictionary is also denoted as {s1, · · · , sK}.
4. K-KSVD method: This method was originally developed in [Nguyen et al., 2013]. It utilizes the kernelized KSVD method to learn the dictionary. In this case, the obtained dictionary is denoted as Φ(Y)A ∈ Rd×K . The default parameter suggested in [Nguyen et al., 2013] are used for this implementation.

As to the proposed method, we fix the regularization parameters λ1 = 0.001 and λ2 = 0.0001 and the maximum iteration number is 100. Please note that the dictionaries in K-KSVD and the proposed method are learned in the feature space, while randDict, K-Means and GRAD methods learn the dictionaries in the input space. Nevertheless, the relation {si}Ki=1 ⊂ {yi}Ni=1 holds for randDict only.

Because the dictionary learning is essentially non-convex problem, for GRAD, K-KSVD and the proposed one, we randomly initialize the dictionary ten times and pick the one with minimum reconstruction error over the training set. In all our experiments, we fix the size of the dictionary to be K = 30. Since all of the compared methods except the proposed one does not equip outlier rejection capability, we resort the classification task to evaluate their performance. We use similar approaches as in [Nguyen et al., 2013] for classification. The algorithm starts with learning a dictionary for each class. In particular, we aggregate all training images for each category into Yc = {yc,1, · · · , yc,Nc }, where c represents the c-th class and Nc is the number of the samples in the c-th class. Let Sc = {sc,1, · · · , sc,K} denotes the learned dictionary in the original space for the c-th class. Such a dictionary can be obtained by using the methods of randDict, KMeans, and GRAD. For K-KSVD and the proposed one, the dictionaries are learned in the feature space and therefore we use Dc = Φ(Yc)Ac to denote the learned kernel dictionary for the c-th class, where Ac is the learned coefficient matrix. Given a query sample z, we first perform the spare coding for each Dc to get the sparse code cc. The sparse setting is the same as the training phase. The reconstruction error for the c-th class is denoted as rc = ||Φ(z) - Φ(Sc)cc||22 for randDict, K-Means and Grad, and rc = ||Φ(z) - Φ(Yc)Accc||22 for K-KSVD and the proposed one. Finally, The test sample is simply classified to the class that gives the smallest reconstruction error.

The classification accuracies on the 4 datasets are reported in Fig.2. It is not surprising that all the curves are monotonically decreasing along the increasing of ρ. From those results we summarize the following observations: (1)K-Means and GRAD do not show significant advantages over randDict. In some especial cases when ρ is large, their performances are even worse than randDict. The main reason is that the dictionary atoms in K-Means and GRAD are learned but not selected and all of atoms are contaminated by the outliers. This means the learning procedures in k-means clustering and gradient learning are strongly affected by the outliers. On the contrary, randDict selects some existing training samples to be the dictionary atoms and therefore some atoms may be clean samples. (2) The performance of the proposed method and K-KSVD is consistently better than that of the other methods. The reason is partially due to the fact that both methods learn dictionaries in the feature space, but not the input space. (3) When the outlier level ρ is small, the performance difference between the proposed method and KKSVD is not significant. However, with the increasing of ρ, the advantage of the proposed method becomes obvious. The intrinsic reason is that the proposed dictionary learning method explicitly isolates the outliers in the training samples. In this sense, the proposed method learns a cleaner dictionary.

Empirically the proposed algorithm works well when the parameter λ2 is in the interval [10-5 , 10-3 ]. When λ2 is too large, the outlier cannot be isolated. On the contrary, when λ2 is too small, too many training samples may be mistakenly isolated as outliers. To study the influence of λ2, we fix λ1 = 0.001 and vary λ2 from 10-6 to 1 and run the algorithm on the 4 datasets. In Fig.3 we show the classification accuracies corresponding to different noise levels. The results show that the proposed algorithm is not very sensitive to the parameter λ2 and a properly-designed robust term indeed plays important role in dictionary learning. In addition, the role of outlier isolation diminishes when λ2 is larger than 10-2 .

5 Conclusions

In this paper, we propose an optimization model for robust kernel dictionary learning which has capability to isolate the outliers. A whole sequence convergent algorithm is developed to solve the non-convex optimization problem and the experimental results show that the proposed robust kernel dictionary learning method provides significant performance improvement.