A New Low-Rank Structurally Incoherent Algorithm for Robust Image Feature Extraction

Abstract

In order to solve the problem in which structurally incoherent low-rank non-negative matrix decomposition (SILR-NMF) algorithms only consider the non-negativity of the data and do not consider the manifold distribution of high-dimensional space data, a new structurally incoherent low rank two-dimensional local discriminant graph embedding (SILR-2DLDGE) is proposed in this paper. The algorithm consists of the following three parts. Firstly, it is vital to keep the intrinsic relationship between data points. By the token, we introduced the graph embedding (GE) framework to preserve locality information. Secondly, the algorithm alleviates the impact of noise and corruption uses the L1 norm as a constraint by low-rank learning. Finally, the algorithm improves the discriminant ability by encrypting the structurally incoherent parts of the data. In the meantime, we capture the theoretical basis of the algorithm and analyze the computational cost and convergence. The experimental results and discussions on several image databases show that the proposed algorithm is more effective than the SILR-NMF algorithm.

1. Introduction

In the past few decades, “curse of dimensionality” [1] problems have been a hot topic for many researchers. A popular dimensionality reduction method, 2D Local Preserving Projections (2DLPP) [2,3,4], has been widely utilized in pattern recognition [5,6,7] and computer vision [8,9,10]. However, on the basis of the 2DLPP unsupervised algorithm, some improved 2DLPP supervised algorithms have also been proposed [11,12,13].

Unfortunately, these above algorithms are similar in that they all use the L2 norm, which will hinder the algorithms’ performance when dealing with noise or abnormal data. Therefore, in order to solve the problem in which the L2 norm is affected by outliers, many researchers use L1-norm-based algorithms as dimension reduction algorithms of distance criteria, which are widely considered to be effective methods [14,15,16,17,18]. For example, L1-PCA [14] and PCA-L1 [15], based on the PCA algorithm, solve the noise and outlier sensitivity in the data through an optimization problem. Finally, the Rotation Invariant L1-norm PCA (R1-PCA) [16] algorithm proposed on the basis of L1-PCA has some PCA properties, which has some PCA properties. In order to solve the general L1 norm local preservation problem, Ref. [17] proposed the LPP algorithm (LPP-L1) based on the PCA-L1 algorithm to more effectively preserve the spatial topological structure. At the same time, in order to solve the outliers and corrosion problems in LPP-L1, a 2DLPP algorithm based on LPP-L1 algorithm (2DLPP-L1) was proposed [18].

In recent years, compared with algorithms based on L1 norm, many feature extraction algorithms using low rank representation (LRR) are conducive to extracting clean information, which has proved to be crucially important when there is noise in the data. The robustness of these algorithms has attracted considerable attention from researchers [19,20,21,22,23,24]. For example, in order to better maintain the lowest rank representation and global structure of data, the LRR in reference [19,20] introduced the single subspace clustering problem into multiple subspace clustering. The subspace structure was recovered from the data damaged by noise or occlusion. On the basis of LRR, a robust PCA (RPCA) was proposed by introducing kernel norm [21,22]. On the basis of LRR, the manifold structure was introduced as the regularization term, and Laplace regularization LRR [23] was proposed and applied to clustering data. Based on the combination of the NNLRS-graph and LRR, the non-negative concept was introduced to propose a non-negative low-rank sparse graph (NNLRS) [24] to maintain the global structure in the data.

A new matrix factorization method with structural incoherence (SILR-NMF) [25] has been introduced into face recognition to improve the recognition performance in the presence of corrupted data. In Ref. [26], a new structurally incoherent method was proposed to improve the additional discrimination ability in face recognition by introducing the structural incoherence. Motivated by [25] and [26], we propose a new algorithm dubbed structural incoherent low rank two-dimensional local discriminant mapping embedding (SILR-2DLDGE), which is superior to the SILR-NMF algorithm who has important limitations. Since the SILR-NMF algorithm can not make full use of the neighborhood relationship between data points. SILR-2DLDGE is implemented in three steps. Firstly, the intra-class weighted matrix graph and the inter-class weighted matrix graph are constructed to maintain the discriminant information of local neighborhood. Secondly, low-rank learning is used to eliminate noise and damage in the data. Finally, the structural incoherencies [25,26] are combined to optimize the discriminant information of different classes.

The four main contributions of this paper are as follows:

• We present a novel algorithm-based 2DLPP, which can simultaneously perform structurally incoherent, optimal graph Laplacian, and low-rank functions in a unified strategy. The algorithm has stronger discriminant ability than SILR-NMF, which fully reveals the structure information of the neighborhood to improve the discrimination ability.

• We propose to introduce intraclass and interclass graphs into the structurally incoherent model to make the data points in the same class more compact and different class data points as far away as possible.

• We used the low-rank feature to ensure that the given data are treated as two parts composed of the low-rank matrix and sparse matrix, representing useful features and nasty noise, respectively, so as to improve the performance of the algorithm using kernel norms as a measure of the regularization term.

• We designed a practical and simple algorithm to cater to the optimization process, and verified the algorithm on five datasets. Actually, our proposed method can achieve better performance.

The rest of study was planned as follows: Section 2 mainly introduces the early algorithms such as LRR, SILR-NMF, and LRMD-SI. Section 3 introduces the details of proposed algorithm, and analyzes the computational complexity and convergence of the SILR-2DLDGE algorithm. In Section 4, we compare the proposed algorithm with other state-of-the-art algorithms, and verify the promising performance of the SILR-2DLDGE algorithm on five databases. The summary and future work can be seen in the last section of the paper.

2. Related Works

2.1. Notation

To facilitate the understanding of this article, we will introduce the related work on some early algorithms, i.e., LRR, LRMD-SI, and SILR-NMF algorithms.

Firstly, we consider $N$ original space samples $X_i\in R^{m\times n}$ into feature space $Y_i\in R^{n\times d}$. Then, the linear transformation shown in Equation (1) can be obtained:

(1)$Y_i=P^T{X}_i,i\in \left(1,2,\cdots ,N\right)$

where $P\in R^{d\times n}$ is a transformation matrix, and $n>d$.

Secondly, we define the matrix $Z{=[\mathrm{z}}_{ij}] \in R^{m\times n}$, and then $Z_i$ or $Z_{.j}$ is used to represent the $i\mathrm{th}$ or $j\mathrm{th}$ row of $Z$, respectively.

Finally, we define the matrix ${\Vert Z\Vert }_1={\displaystyle {\sum }_{i,j}\left|z_{ij}\right|}$, ${\Vert Z\Vert }_F=\sqrt{{\displaystyle {\sum }_{}{\Vert Z_{i.}\Vert }_2^2}}$, and ${\Vert Z\Vert }_{2,1}={{\displaystyle {\sum }_j\Vert Z_{.j}\Vert }}_2$, respectively. Additionally, ${\Vert Z\Vert }_{*}$ represents the kernel norm which is used to compute the sum of the singular matrix.

2.2. LRR

The LRR [19] algorithm is different from sparse representation [27]. LRR adopts a joint approach to transfer the recovery of damaged data to multiple subspaces, which is an effective subspace segmentation algorithm. At the same time, the global structure of data is obtained by searching the lowest rank representation of all data.

Each data vector in $X$ can be learned by the linear combination of the $T$:

(2)$X=TZ,$

where $Z=\left[Z_1,Z_2,\dots Z_n\right]$ is the coefficient matrix being, representation $X$:

(3)$\underset{Z}{\min }{\Vert Z\Vert }_{*}^{}, \mathrm{s}.\mathrm{t}. \mathrm{X}=TZ,$

Here, ${\Vert \cdot \Vert }_{*}$ denotes the nuclear norm.

The noise E can be separated from corrupted data as follows [19]:

(4)$\underset{Z}{\min }{\Vert Z\Vert }_{*}^{}+\lambda {\Vert E\Vert }_{ 2,1} \mathrm{s}.\mathrm{t}. Z=\mathrm{X}-TE,$

where $\lambda$ is a parameter.

2.3. SILR-NMF

In view of the fact that face is often affected by noise, SILR-NMF [25] is a better choice for feature representation and image classification. We can obtain Equation (5) as the final form of its objective function:

(5)$\underset{Z,E,U\ge 0,V\ge 0}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert X-UV\Vert }_F^2+}{\displaystyle \sum _{i\ne j} \alpha {\Vert Z_i\Vert }_{*}+\beta {\Vert E_i\Vert }_1}+\gamma {\displaystyle \sum _{i\ne j}{\Vert Z_j^T{Z}_i\Vert }_F^2} s.t. Z=X-E$

where $\alpha$, $\beta$, and $\gamma$ are positive parameters.

2.4. LRMD-SI

In view of the fact that face is often affected by noise, Wei et al. [26] introduced a structurally incoherent constraint based on low-rank matrix decomposition (LRMD-SI), which is considered to have good performance in image representation and classification. We can obtain Equation (6) as the final form of its objective function:

(6)$\underset{Z,E}{{\displaystyle \min }} {\displaystyle \sum _{i=1}^c\left\{{\Vert Z_i\Vert }_{*}+\lambda {\Vert E_i\Vert }_1\right\}+\eta {\displaystyle \sum _{i\ne j}{\Vert Z_j^T{Z}_i\Vert }_F^2}} s.t. X_i=E_i+Z_i,$

where matrix ${\Vert .\Vert }_F^{}$ is the F-norm of the matrix, $\eta$ is the parameter, and $i\in \left(1,2,\dots ,c.\right)$.

3. Structurally Incoherent Low-Rank 2D Local Discriminant Graph Embedding

This section discusses the objective function of the SILR-2DLDGE algorithm, and describes its optimization process, computational complexity, and convergence.

3.1. The Objective Function of SILR-2DLDGE

In low-rank matrix recovery, the training dataset $X=\left[{\tilde{X}}_1,{\tilde{X}}_2,\cdots ,{\tilde{X}}_C\right]$, which is divided into the low-rank matrix $Z=\left[{\tilde{Z}}_1,{\tilde{Z}}_2,\cdots ,{\tilde{Z}}_C\right]$ and the noise matrix $E=\left[{\tilde{E}}_1,{\tilde{E}}_2,\cdots ,{\tilde{E}}_C\right]$, where ${\tilde{X}}_i=\left[{\tilde{X}}_1,{\tilde{X}}_2,\cdots ,{\tilde{X}}_N\right]$$\left(i=1,2,\cdots ,C\right)$, $C$ is the number of class of training sample $X$, and N is the total number of data in the $i$ class of training sample $X$. Then, we have the following:

(7)$\underset{{\tilde{Z}}_i,{\tilde{E}}_i}{{\displaystyle \min }} {\displaystyle \sum _{i=1}^C\left\{{\Vert {\tilde{Z}}_i\Vert }_{*}+\lambda {\Vert {\tilde{E}}_i\Vert }_1\right\}} s.t. {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i$

where $\lambda$ is a parameter.

The following equation can be obtained by encoding the 2DLDGE algorithm into Equation (10) and introducing an orthogonal constraint on $Y$:

(8)$$\begin{gathered} \underset{{\tilde{Z}}_i,{\tilde{E}}_i,{\tilde{P}}_i}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}{\displaystyle \sum _{i=1}^C\left\{\alpha {\Vert {\tilde{Z}}_i\Vert }_{*}+\beta {\Vert {\tilde{E}}_i\Vert }_1\right\}} \\ s.t. {\tilde{Z}}_i={\tilde{B}}_i, {\tilde{Y}}_i={\tilde{P}}_i{}^T{\tilde{Z}}_i and {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i \end{gathered}$$

where the balance parameters $\alpha >0$ and $\beta >0$.

We define the scatter matrix $S_{ij}^w$ (intra-class) and the scatter matrix $S_{ij}^b$ (inter-class) in Equation (8):

(9)$S_{ij}^w=\{\begin{array}{l}1,\hspace{0.17em}if \hspace{0.17em}{X}_i\in N_{K_c}^{+}(X_j)\hspace{0.17em}or\hspace{0.17em}{X}_j\in N_{K_c}^{+}(X_i) \\ 0,\hspace{0.17em}otherwise.\end{array}\hspace{0.17em}$

(10)$S_{ij}^b=\{\begin{array}{l}1,\hspace{0.17em}{X}_i\in N_{K_b}^{+}(X_j)\hspace{0.17em}or\hspace{0.17em}{X}_j\in N_{K_b}^{+}(X_i)\\ 0,\hspace{0.17em}otherwise.\end{array}\hspace{0.17em}$

where $N_{K_c}^{+}(X_i)$ or $N_{K_c}^{+}(X_j)$ in the same class means the index set of the K_c nearest neighbors of the sample $X_i$ or $X_j$, respectively. Additionally, $N_{K_b}^{+}(X_i)$ or $N_{K_b}^{+}(X_j)$ in the different class express the index set of the K_b nearest neighbors of the sample $X_i$ or $X_j$, respectively.

We added structural incoherence to Equation (8), which is expressed as follows to improve the discrimination ability and further distinguish the different classes of the Equation (7) algorithm:

(11)$$\begin{gathered} \underset{{\tilde{Z}}_i,{\tilde{E}}_i,{\tilde{P}}_i}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+} {\displaystyle \sum _{i=1}^C\left\{\alpha {\Vert {\tilde{Z}}_i\Vert }_{*}+\beta {\Vert {\tilde{E}}_i\Vert }_1\right\}}+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{Z}}_i\Vert }_F^2} \\ s.t. {\tilde{Z}}_i={\tilde{B}}_i, {\tilde{Y}}_i={\tilde{P}}_i{}^T{\tilde{Z}}_i and {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i \end{gathered}$$

In Equation (11), the algorithm improves the incoherence of algorithm structure by adding the F-norm between different pairs of ${\tilde{Z}}_i$ and ${\tilde{Z}}_j$. The first term is to maintain the neighborhood information in the low dimensional subspace consistent with the clean data structure. The second term is used to learn a low rank matrix to ensure that its noise interference is reduced as much as possible. The last term has the discrimination ability and further distinguishes the different classes in the data.

3.2. The Optimization of SILR-2DLDGE

Due to the product term ${\tilde{Z}}_i^T$${\tilde{Z}}_j$, the minimization Equation (11) is non-convex. We iteratively replaced low-rank matrices ${\tilde{Z}}_i$ to solve the following problems:

(12)$$\begin{gathered} \underset{{\tilde{Z}}_i,{\tilde{E}}_i,{\tilde{P}}_i}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}\alpha {\Vert {\tilde{Z}}_i\Vert }_{*}+ \beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{Z}}_i\Vert }_F^2} \\ s.t. {\tilde{Z}}_i={\tilde{B}}_i, {\tilde{Y}}_i={\tilde{P}}_i{}^T{\tilde{Z}}_i and {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i \end{gathered}$$

By introducing the auxiliary variable ${\tilde{B}}_i$ to solve Equation (12), the following equation can be obtained:

(13)$$\begin{gathered} \underset{{\tilde{Z}}_i,{\tilde{E}}_i,{\tilde{P}}_i}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}\alpha {\Vert {\tilde{Z}}_i\Vert }_{*}+\beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2} \\ s.t. {\tilde{Z}}_i={\tilde{B}}_i, {\tilde{Y}}_i={\tilde{P}}_i{}^T{\tilde{Z}}_i and {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i \end{gathered}$$

Augmented Lagrange multipliers or ADMM [28] and other techniques can be used to solve Equation (13). We solved Equation (13) by performing singular value decomposition (SVD) at each iteration step, but it may have high computational complexity. We used the nuclear norm property [29] to reduce the complexity of Equation (13) as follows [30]:

(14)${\Vert {\tilde{Z}}_i\Vert }_{*}=\underset{R.H}{\min }\frac{1}{2}\left({\Vert H\Vert }_F^2+{\Vert R\Vert }_F^2\right), s.t. {\tilde{Z}}_i=RH$

where ${\tilde{Z}}_i\in R^{a\times bN}$ is an arbitrary matrix. If $r\ge \min \left(a,bN\right)$ of ${\tilde{A}}_i$, the factorization ${\tilde{Z}}_i=R_{a\times r}{H}_{r\times bN}$ is obtained by the minimum solution of the above equation.

We can write the equivalent of Equation (13) using the conclusion of Equation (14) as follows:

(15)$$\begin{gathered} \underset{{\tilde{Z}}_i,{\tilde{E}}_i,\tilde{P}}{{\displaystyle \min }} {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}\frac{\alpha }{2}\left({\Vert R\Vert }_F^2+{\Vert H\Vert }_F^2\right)+\beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2} \\ s.t. {\tilde{Z}}_i={\tilde{B}}_i, {\tilde{Y}}_i={\tilde{P}}_i{}^T{\tilde{Z}}_i and {\tilde{X}}_i={\tilde{Z}}_i+{\tilde{E}}_i \end{gathered}$$

The augmented Lagrange function of Equation (15) is detailed as follows:

(16)$\begin{array}{l}L\left({\tilde{Z}}_i,{\tilde{B}}_i,\mu ,{\tilde{E}}_i,{\tilde{P}}_i,R,H,M_1,M_2,M_3\right)\\ = {\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}\frac{\alpha }{2}\left({\Vert R\Vert }_F^2+{\Vert H\Vert }_F^2\right)+\beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2} +tr\left(M_1^T\left({\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i\right)\right)\\ +tr\left(M_2^T\left({\tilde{Z}}_i-{\tilde{B}}_i\right)\right)+tr\left(M_3^T\left({\tilde{Z}}_i-RH\right)\right) +\frac{\mu }{2}\left({\Vert {\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i\Vert }_F^2+{\Vert {\tilde{Z}}_i-{\tilde{B}}_i\Vert }_F^2+{\Vert {\tilde{Z}}_i-RH\Vert }_F^2\right)\\ ={\displaystyle \sum _{i,j}{\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^w-\gamma {\Vert Y_i-Y_j\Vert }_2^2{S}_{ij}^b+}\frac{\alpha }{2}\left({\Vert R\Vert }_F^2+{\Vert H\Vert }_F^2\right)+\beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2}+\frac{\mu }{2}{\Vert {\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i+\frac{M_1}{\mu }\Vert }_F^2\\ +\frac{\mu }{2}{\Vert {\tilde{Z}}_i-{\tilde{B}}_i+\frac{M_2}{\mu }\Vert }_F^2 +\frac{\mu }{2}{\Vert {\tilde{Z}}_i-RH+\frac{M_3}{\mu }\Vert }_F^2-\frac{1}{2\mu }\left({\Vert M_1\Vert }_F^2+{\Vert M_2\Vert }_F^2+{\Vert M_3\Vert }_F^2\right)\\ =2{\tilde{P}}_i{}^T{\tilde{Z}}_i\left(L^w\otimes I\right){\tilde{Z}}_i{}^T{\tilde{P}}_i-2{\tilde{P}}_i{}^T{\tilde{Z}}_i\left(L^b\otimes I\right){\tilde{Z}}_i{}^T{\tilde{P}}_i+\frac{\alpha }{2}\left({\Vert R\Vert }_F^2+{\Vert H\Vert }_F^2\right)+\beta {\Vert {\tilde{E}}_i\Vert }_1+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2}\\ +\frac{\mu }{2}{\Vert {\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i+\frac{M_1}{\mu }\Vert }_F^2 +\frac{\mu }{2}{\Vert {\tilde{Z}}_i-{\tilde{B}}_i+\frac{M_2}{\mu }\Vert }_F^2 +\frac{\mu }{2}{\Vert {\tilde{Z}}_i-RH+\frac{M_3}{\mu }\Vert }_F^2-\frac{1}{2\mu }\left({\Vert M_1\Vert }_F^2+{\Vert M_2\Vert }_F^2+{\Vert M_3\Vert }_F^2\right)\end{array}$

where the penalty parameter is $\mu >0$, and the Lagrange multipliers are $M_1$, $M_2$, and $M_3$.

Next, we separated the solutions for each variable.

• (1). Fix ${\tilde{Z}}_i$, ${\tilde{E}}_i$, $R$, $H$, and $P$, and update ${\tilde{B}}_i$: since ${\tilde{Z}}_i$, ${\tilde{E}}_i$, $R$, $H$, and $P$ are fixed, we can obtain the following equation from Equation (16):

(17) $\frac{\mu }{2}{\Vert {\tilde{Z}}_i-{\tilde{B}}_i+\frac{M_2}{\mu }\Vert }_F^2+\eta {\displaystyle \sum _{i\ne j}{\Vert {\tilde{Z}}_j^T{\tilde{B}}_i\Vert }_F^2}$

Letting the partial derivative of variable B be zero in Equation (17):

(18)$\mu \left({\tilde{B}}_i-{\tilde{Z}}_i-\frac{M_2}{\mu }\right)+2\eta {\displaystyle \sum _{i\ne j}{\tilde{Z}}_j^{}{\tilde{Z}}_j^T{\tilde{B}}_i}=0$

Then, we have the following:

(19)${\tilde{B}}_i={\left(2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_j^T}+\mu I\right)}^{-1}\left(\mu {\tilde{Z}}_i+M_2\right)$

• (2). Fix ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$, $H$, and $P$, and update ${\tilde{Z}}_i$: by fixing other variables, we can obtain the following equation from Equation (16):

(20) $\begin{array}{l}2{\tilde{P}}_i{}^T{\tilde{Z}}_i\left(L^w\otimes I\right){\tilde{Z}}_i{}^T{\tilde{P}}_i-2{\tilde{P}}_i{}^T{\tilde{Z}}_i\left(L^b\otimes I\right){\tilde{Z}}_i{}^T{\tilde{P}}_i+\\ {\Vert {\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i+\frac{M_1}{\mu }\Vert }_F^2+{\Vert {\tilde{Z}}_i-{\tilde{B}}_i+\frac{M_2}{\mu }\Vert }_F^2+{\Vert {\tilde{Z}}_i-RH+\frac{M_3}{\mu }\Vert }_F^2\end{array}$

In Equation (20), letting the partial derivative of ${\tilde{Z}}_i$ be zero, we have:

(21)$4{\tilde{P}}_i{\tilde{P}}_i{}^T{\tilde{Z}}_i\left(\left(L^w-L^b\right)\otimes I\right)+3{\tilde{Z}}_i-\left({\tilde{Q}}_1+{\tilde{Q}}_2+{\tilde{Q}}_3\right)=0$

where ${\tilde{Q}}_3=RH-\frac{M_3}{\mu }$, ${\tilde{Q}}_2={\tilde{B}}_i-\frac{M_2}{\mu }$ and ${\tilde{Q}}_1={\tilde{X}}_i-{\tilde{E}}_i+\frac{M_1}{\mu }$.

• (3). Fix ${\tilde{Z}}_i$, ${\tilde{B}}_i$, $R$, $H$, and $P$, and update ${\tilde{E}}_i$: to fix other variables, the optimal ${\tilde{E}}_i$ can be obtained in Equation (16); thus, we have:

(22) $\beta {\Vert {\tilde{E}}_i\Vert }_1+\frac{\mu }{2}{\Vert {\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i+\frac{M_1}{\mu }\Vert }_F^2$

To simplify the solution of the above problem, we introduce a shrinkage operator according to reference [31]. In addition, $S_{\epsilon }\left[X\right]=\max \left(\left|X\right|-\epsilon ,0\right).\mathrm{sgn}(X)$ represents the soft threshold to achieve the above purposes [32]. Finally, it is converted into Equation (23):

(23)${\tilde{E}}_i=S_{\frac{\beta }{\mu }}\left({\tilde{X}}_i-{\tilde{Z}}_i+\frac{M_1}{\mu }\right)$

• (4). Fix ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $H$, and $P$, and update $R$: we obtain the following Equation (24) by fixing all variables except $R$ in Equation (16):

(24) $\frac{1}{2}\left(\alpha {\Vert R\Vert }_F^2+\mu {\Vert {\tilde{A}}_i-RH+\frac{M_3}{\mu }\Vert }_F^2\right)$

In Equation (24), letting the partial derivative of variable $R$ be zero, we obtain:

(25)$\begin{array}{l}\mu RH{H}^T+\alpha R-\left(M_3+\mu {\tilde{Z}}_i\right)H^T=0\\ \Rightarrow R=\left(M_3+\mu {\tilde{Z}}_i\right)H{\left(\mu RH{H}^T+\alpha R\right)}^{-1}\end{array}$

• (5). Fix ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$, and $P$, and update $H$: we fix all other variables in Equation (16) except $H$; thus, the following equation can be obtained:

(26) $\frac{1}{2}\left(\alpha {\Vert H\Vert }_F^2+\mu {\Vert {\tilde{Z}}_i-RH+\frac{M_3}{\mu }\Vert }_F^2\right)$

The following equation is obtained to use the partial derivative of Equation (26):

(27)$\begin{array}{l}\mu R^{T}RH+\alpha H-\left(\mu {\tilde{Z}}_i+M_3\right)H^T=0\\ \Rightarrow H={\left(\alpha I+\mu R^{T}R\right)}^{-1}{R}^T\left(\mu {\tilde{Z}}_i+M_3\right)\end{array}$

• (6). Fix ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$, and $H$, and update $P$: we simplify Equation (16) and obtain the solution of variable $P$:

(28) $\mathrm{\Psi } =\arg \underset{{\tilde{P}}_i}{\min } 2{\tilde{P}}_i{}^T\tilde{Z}\left(L^w\otimes I\right){\tilde{Z}}^T{\tilde{P}}_i-2{\tilde{P}}_i{}^T\tilde{Z}\left(L^b\otimes I\right){\tilde{Z}}^T{\tilde{P}}_i$

(29) $$\begin{gathered} \mathrm{Equation} (28) \mathrm{will} \mathrm{be} \mathrm{converted} \mathrm{as} \mathrm{follows}: \\ 2{\tilde{P}}_i{}^T\tilde{Z}\left(\left(L^w-L^b\right)\otimes I\right){\tilde{Z}}^T{\tilde{P}}_i= 2{\tilde{P}}_i{}^T(\tilde{X}-\tilde{E})\left(\left(L^w-L^b\right)\otimes I\right){(\tilde{X}-\tilde{E})}^T{\tilde{P}}_i \end{gathered}$$

In order to solve the equation better, we add a constraint adjustment; the equation is as follows:

(30)${\tilde{P}}_i{}^T\left(D^w-D^b\right){(\tilde{X}-\tilde{E})}^T{\tilde{P}}_i= 1$

Equations (29) and (30) can be transformed into the following formula:

(31)$\underset{p^T\left(X-E\right)\left(D^w-D^b\right){\left(X-E\right)}^{T}p=1}{\arg \min } {\tilde{P}}^T\left(X-E\right)\left(\left(L^w-L^b\right)\otimes I\right){\left(X-E\right)}^T\tilde{P}$

The solution of Equation (31) is the same as solving the generalized eigenvalue equation as follows:

(32)$\begin{array}{l}\left(X-E\right)\left(\left(L^w-L^b\right)\otimes I\right){\left(X-E\right)}^{T}p\\ =\lambda \left(X-E\right)\left(D^w-D^b\right){\left(X-E\right)}^{T}p\end{array}$

Algorithm 1 gives the concrete steps of SILR-2DLDGE.

Algorithm 1 SILR-2DLDGE

Input: Parameters $\alpha$, $\beta$, $\gamma$ and $\eta$; and Training set X in Equation (14).Initialization: ${\tilde{E}}_i$ = 0, ${\tilde{B}}_i$ = 0, $\rho$ > 0, $M_1$ = 0, $M_2$ = 0, $M_3$ = 0. RepeatStep 1. Fixing ${\tilde{Z}}_i$, ${\tilde{E}}_i$, $R$, $H$ and $P$, update ${\tilde{B}}_i$ with Equation (17); Step 2. Fixing ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$, $H$ and $P$, update ${\tilde{Z}}_i$ with Equation (20); Step 3. Fixing ${\tilde{Z}}_i$, ${\tilde{B}}_i$, $R$, $H$ and $P$, update ${\tilde{E}}_i$ with Equation (22); Step 4. Fixing ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $H$ and $P$, update $R$ with Equation (24);Step 5. Fixing ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$ and $P$, update $H$ with Equation (26);Step 6. Fixing ${\tilde{Z}}_i$, ${\tilde{B}}_i$, ${\tilde{E}}_i$, $R$ and $H$, update $P$ with Equation (28); Step 7. Update parameter $M_1$, $M_2$, and $M_3$ as follows:Update $M_1$ through $M_1=\mu \left({\tilde{X}}_i-{\tilde{A}}_i-{\tilde{E}}_i\right)+M_1$; Update $M_2$ through $M_2=\mu \left({\tilde{A}}_i-{\tilde{B}}_i\right)+M_2$; Update $M_3$ through $M_3=\mu \left({\tilde{A}}_i-RH\right)+M_3$; Step 8. Update $\mu$ through $\mu =\min (\rho \mu ,\max \mu )$; Step 9. Update $t$ through t = t + 1; Until Equation (11) is converged.Step 10. Obtain the (${\tilde{B}}_i$, ${\tilde{A}}_i$, ${\tilde{E}}_i$, $R$, $H$, $P$) optimal solution;Output: Obtain projection matrix $P$

The flow chart of the SILR-2DLDGE method is shown in Figure 1.

3.3. The Convergence Analysis

Firstly, we analyzed and proved the weak convergence of the proposed SILR-2DLDGE algorithm. In some cases, any limit point of the iterative sequence generated by the SILR-2DLDGE algorithm is a stationary point satisfying the Karush–Kuhn–Tucker (KKT) condition [33].

Assuming that the SILR-2DLDGE algorithm reaches a stationary point, any convergence point must satisfy the KKT condition, which is a necessary condition for a local optimal solution. Therefore, we can derive the KKT condition of Equation (11) as follows (note that Lagrange multiplier is not involved in the process of solving $P$; thus, we do not prove its KKT condition):

(33)$$\begin{gathered} {\tilde{Z}}_i-RH=0, {\tilde{Z}}_i-{\tilde{B}}_i=0,{\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i=0 \\ \frac{\partial L}{\partial {\tilde{Z}}_i}= 4p{p}^T{\tilde{Z}}_i\left(\left(L^w-L^b\right)\otimes I\right)+3{\tilde{Z}}_i-\left({\tilde{Q}}_1+{\tilde{Q}}_2+{\tilde{Q}}_3\right)=0 \\ \frac{\partial L}{\partial {\tilde{B}}_i}= \mu \left({\tilde{B}}_i-{\tilde{Z}}_i-\frac{M_2}{\mu }\right)+2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_j^T{\tilde{B}}_i}=0 \\ \frac{\partial L}{\partial R}=\alpha R-\left(\mu {\tilde{Z}}_i+M_3\right)H^T+\mu RH{H}^T=0 \\ \begin{array}{l}\left(X-E\right)\left(\left(L^w-L^b\right)\otimes I\right){\left(X-E\right)}^{T}p\\ =\lambda \left(X-E\right)\left(D^w-D^b\right){\left(X-E\right)}^{T}p\end{array} \\ \frac{\partial L}{\partial H}=\alpha H-\left(\mu {\tilde{Z}}_i+M_3\right)H^T+\mu R^{T}RH=0 \\ {M}_1\in \beta {\partial }_{{\tilde{E}}_i}{\Vert {\tilde{E}}_i\Vert }_1 \end{gathered}$$

It can be deduced from the last relationship in Equation (33) that:

(34)$\frac{M_1}{\mu }+{\tilde{X}}_i-{\tilde{Z}}_i\in \beta \frac{{\partial }_{{\tilde{E}}_i }{\Vert {\tilde{X}}_i-{\tilde{Z}}_i\Vert }_1}{\mu }+{\tilde{X}}_i-{\tilde{Z}}_i\triangleq {{\displaystyle \mathrm{\Lambda }}}_{\frac{\beta }{\mu }}\left({\tilde{X}}_i-{\tilde{A}}_i\right)$

where the application element to ${\tilde{X}}_i-{\tilde{Z}}_i$ is the scalar function ${{\displaystyle \mathrm{\Lambda }}}_{\frac{\beta }{\mu }}\left(t\right)=t+\frac{\beta }{\mu }\partial \left|t\right|$. From [34], the following relationships will be learned:

(35)${\tilde{E}}_i={\mathrm{\Lambda }}_{\frac{\beta }{\mu }}^{-1}\left(\frac{M_1}{\mu }+{\tilde{X}}_i-{\tilde{Z}}_i\right)=S\left(\frac{M_1}{\mu }+{\tilde{X}}_i-{\tilde{Z}}_i,\frac{\beta }{\mu }\right)$

where $S\left(x,\epsilon \right)=\max \left(\left|x\right|-\epsilon ,0\right).sign\left(x\right)$. Thus, the following equation can be obtained with the KKT condition:

(36)$$\begin{gathered} {\tilde{Z}}_i-RH=0, {\tilde{Z}}_i-{\tilde{B}}_i=0,{\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i=0 \\ \frac{\partial L}{\partial {\tilde{Z}}_i}= 4p{p}^T{\tilde{Z}}_i\left(\left(L^w-L^b\right)\otimes I\right)+3{\tilde{Z}}_i-\left({\tilde{Q}}_1+{\tilde{Q}}_2+{\tilde{Q}}_3\right)=0 \\ \frac{\partial L}{\partial {\tilde{B}}_i}= \mu \left({\tilde{B}}_i-{\tilde{Z}}_i-\frac{M_2}{\mu }\right)+2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_j^T{\tilde{B}}_i}=0 \\ \frac{\partial L}{\partial R}=\alpha R+\mu RH{H}^T-\left(\mu {\tilde{Z}}_i+M_3\right)H^T=0 \\ \frac{\partial L}{\partial H}=\alpha H+\mu R^{T}RH-\left(\mu {\tilde{Z}}_i+M_3\right)H^T=0 \\ {\tilde{E}}_i=S\left({\tilde{X}}_i-{\tilde{Z}}_i+\frac{M_1}{\mu },\frac{\beta }{\mu }\right) \end{gathered}$$

On the basis of the above conditions, the point convergence satisfying the KKT condition is proved.

Theorem 1 [25].

Suppose that ${\left\{{\theta }^j\right\}}_{j=1}^{\infty }$is bounded to be generated by the SILR-2DLDGE algorithm and let$\theta$$\triangleq$$\left({\tilde{Z}}_i,{\tilde{B}}_i,H,R,{\tilde{E}}_i,{\tilde{P}}_i,M_1,M_2,M_3\right)$. If${\lim }_{j\to \infty }\left\{{\theta }^{j+1}-{\theta }^j\right\}=0$satisfies KKT conditions and converges, then it can be obtained that any point converges to the KKT point in${\left\{{\theta }^j\right\}}_{j=1}^{\infty }$.

Proof of Theorem 1.

$M_i^{+}$ is the next point of $M_i^{}$ in the sequence ${\left\{M_i^j\right\}}_{j=1}^{\infty }$. Firstly, $M_1$, $M_2$, and $M_3$ can be obtained, as shown in the following equation:

(37) $$\begin{gathered} M_1^{+}\leftarrow \mu \left({\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i\right)+M_1 \\ {M}_2^{+}\leftarrow \mu \left({\tilde{Z}}_i-{\tilde{B}}_i\right)+M_2 \\ {M}_3^{+}\leftarrow \mu \left({\tilde{Z}}_i-RH\right)+M_3 \end{gathered}$$

If $\left({\tilde{X}}_i-{\tilde{Z}}_i-{\tilde{E}}_i\right)\to 0$, $\left({\tilde{Z}}_i-{\tilde{B}}_i\right)\to 0$ and $\left({\tilde{Z}}_i-RH\right)\to 0$, then we can obtain that $\left(M_1^{+}-M_1\right)\to 0$, $\left(M_2^{+}-M_2\right)\to 0$ and $\left(M_3^{+}-M_3\right)\to 0$. Thus, the sequence variables of ${\left\{M_1^j\right\}}_{j=1}^{\infty }$, ${\left\{M_2^j\right\}}_{j=1}^{\infty }$ and ${\left\{M_3^j\right\}}_{j=1}^{\infty }$ converge to the stationary point, and we obtain the first three conditions in Equation (37).

Next, according to the SILR-2DLDGE algorithm and the fourth KKT condition, we obtain Equation (38):

(38)$\begin{array}{l}\left(2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_{\mathit{ij}}^T}+\mu I\right)\left({\tilde{B}}_i^{+}-{\tilde{B}}_i\right)\\ = \left(\mu {\tilde{Z}}_i+M_2\right)-\left(2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_{\mathit{ij}}^T}+\mu I\right){\tilde{B}}_i\end{array}$

When $\left({\tilde{B}}_i^{+}-{\tilde{B}}_i\right)\to 0$, we can derive $M_2-2\eta {\displaystyle \sum _{j\ne i}{\tilde{Z}}_j^{}{\tilde{Z}}_j^T{}_i}{\tilde{B}}_i\to 0$ from the second condition.

The following equation is obtained by the third KKT condition:

(39)$\begin{array}{l}2\mu \left({\tilde{Z}}_i^{+}-{\tilde{Z}}_i\right)\\ =\mu \left(Q_1+Q_2+Q_3\right)-4p{p}^T{\tilde{Z}}_i\left(\left(L^w-L^b\right)\otimes I\right)-2\mu {\tilde{Z}}_i\end{array}$

Thus, $4p{p}^T{\tilde{Z}}_i\left(\left(L^w-L^b\right)\otimes I\right)+2\mu {\tilde{Z}}_i-\mu \left(Q_1+Q_2+Q_3\right)\to 0$ when $\left({\tilde{Z}}_i^{+}-{\tilde{Z}}_i\right)\to 0$.

Similarly, the following equation can be obtained:

(40)$\begin{array}{l}\left(R^{+}-R\right)\left(\mu H{H}^T+\alpha I\right)\\ =\mu \left({\tilde{Z}}_i+\frac{M_2}{\mu }\right)H^T-R\left(\mu H{H}^T+\alpha I\right)\end{array}$

When $\left(R^{+}-R\right)\to 0$, $\alpha R+\mu RH{H}^T-\mu \left({\tilde{Z}}_i+\frac{M_2}{\mu }\right)H^T\to 0$ can be obtained. Similarly:

(41)$\begin{array}{l}\left(\alpha I+\mu R^{T}R\right)\left(H^{+}-H\right)\\ =\mu R^T\left({\tilde{Z}}_i+\frac{M_2}{\mu }\right)-\left(\alpha I+\mu R{R}^T\right)H\end{array}$

When $\left(H^{+}-H\right)\to 0$, we have $\alpha H+\mu R^{T}RH-\mu R^T\left({\tilde{Z}}_i+\frac{M_2}{\mu }\right)\to 0$.

From the last condition in Equation (35), the following equation is as follows:

(42)${\tilde{E}}_i^{+}-{\tilde{E}}_i=S\left({\tilde{X}}_i-{\tilde{Z}}_i+\frac{M_1}{\mu },\frac{\beta }{\mu }\right)-{\tilde{E}}_i$

We can obtain the last condition when ${\tilde{E}}_i^{+}-{\tilde{E}}_i\to 0$.

Both sides of Equations (38)–(42) indicate tending to 0 when ${\lim }_{j\to \infty }\left\{{\theta }^{j+1}-{\theta }^j\right\}=0$. Finally, the sequence variable of ${\left\{{\theta }^j\right\}}_{j=1}^{\infty }$ can be obtained asymptotically to satisfy the KKT condition in Equation (11).

QED □

3.4. Computational Complexity

The main computation steps of SILR-2DLDGE algorithm are shown in steps 2 and 6. In the SILR-2DLDGE algorithm, the complexity of step 2 is O(t(2a²l)) to solve the Sylvester equation problem, and the complexity of step 6 is O(t(2a³)). Therefore, we can obtain the total complexity of the SILR-2DLDGE algorithm, which is O(t(2(a³ + a²l))), where t is the number of iterations.

4. Experiment Results

In this section, we introduce several representative experiments to verify the robustness of the proposed SILR-2DLDGE algorithm. We compared the results of the proposed algorithm with 2DLPP [2,3,4], LPP-L1 [17], 2DLP-L1 [18], LRR [19,20], SILR-NMF [25], LRMD-SI [26], and LPP [35,36] algorithms on FERET [37], ORL [38], coil100 [39], Yale [40], AR [41], and PolyU [42] databases, respectively.

4.1. Selection of Samples and Parameters

In each subsequent experiment, we will randomly select T = 5, 6, 10, 5, 6, 3 samples from each class in FERET, ORL, COIL 100, Yale, AR, and PolyU databases for training. In each running experiment, NN classifier was used for classification, and the experiment was repeated 10 times.

In all the following experiments, we used an iterative algorithm to derive solutions of the L1-norm algorithms and set the maximum value of iterations of the relevant iterative algorithm to 500. The LPP, 2DLPP, LPP-L1, and 2DLPP-L1 algorithms based on local graph embedding are obtained using k-nearest neighbors, where k = T − 1 (T is the training sets of each databases) is well collected in the observation space [43]. We choose the parameters of LRR, LRMD-SI, and SILR-NMF models in the references. From the objective function of the SILR-2DLDGE algorithm in Equation (13), the parameter range of α, β, and η were selected from [0.001, 0.01, 0.1, 1, 10, 100, 1000], and the γ parameter was selected from [0.1, 0.2, ……, 0.9, 1] to evaluate the effect, the values of parameters in our algorithm are given in Table 1.

4.2. Experiments on Occlusion Databases

To verify the robustness of the algorithm, we will randomly add 10 × 10 blocks to different positions of images to carry out continuous occlusion experiments in FERET, ORL, and COIL100 databases.

• FERET database

The FERET database is mainly used to study changes in pose, illumination, and facial expression. There are 200 classes in the FERET database, and each class has seven images with a resolution of 40 × 40 pixels, resulting in a total of 1400 gray-scale images.

• ORL database

The ORL database is mainly used to study changes with different expressions, postures, and illuminations. There are 40 classes in the ORL database., and each class has 10 images with a resolution of 56 × 46 pixels, resulting in a total of 400 gray-scale images.

• COIL 100 database

The COIL 100 object database is mainly used to study changes with different illuminations. There are 100 subjects, and each subject has 72 images with a resolution of 32 × 32 pixels, resulting in a total of 7200 gray-scale images. Figure 2 shows some occlusion images on the three different databases.

To verify the effectiveness of continuous occlusion in FERET, ORL, and COIL100 databases, respectively, we randomly added 10 × 10 blocks to different positions of images. Figure 3, Figure 4 and Figure 5 show the average recognition rates (%) of different dimension changes on the three different databases.

We ran a set of experimental results and compared them with seven methods; the results are presented in this section. According to the settings in Section 4.1, when each class in FERET, ORL, and COIL100 databases randomly selects T sample points to form a training sample set, the optimal average recognition rate (%) of the seven algorithms is shown in Table 2.

4.3. Experiments on Noise Databases

Experiments were carried out on different levels of random pixel corruptions to verify the robustness of the algorithm. “Salt & pepper” noise with a density of 0.1 was added to the Yale, AR, and PolyU databases for corrosion experiments.

• Description of the Yale database

The Yale face database is mainly used to study changes in facial expressions and lighting conditions. There are 15 classes in the Yale database, and each class has 11 images with a resolution of 50 × 40 pixels, resulting in a total of 165 gray-scale images.

• Description of the AR database

The AR database is mainly used to study the changes in lighting conditions, facial expressions, and occlusion. There are 70 men and 56 women with a total of 126 people in the AR database, including 4000 color images with a resolution of 50 × 40 pixels.

• Description of the PolyU palmprint database

The PolyU database is mainly used to study image changes in two periods. There are 100 different palms in the PolyU database, and each palm has six samples, resulting in a total of 600 gray images with a resolution of 64 × 64 pixels. Figure 6 shows the original images and some corrosion images on the three different databases.

For continuous corrosion experiments, we randomly added a density of 0.1 of “salt & pepper” noise to the Yale, AR, and PolyU databases. Then, we ran a set of experiments and compared them with seven methods to evaluate the proposed SILR-2DLDGE algorithm. Figure 7, Figure 8 and Figure 9 show the average recognition rates (%) of different dimension changes on the three different databases.

As in the last experiment, we can obtain the best average recognition rates (%) of the eight algorithms according to the settings in Section 4.1, as shown in Table 3.

4.4. Convergence Study

In the last experiment, we will further prove the convergence of the SILR-2DLDGE algorithm on the FERET and AR databases, respectively. The objective function values decreased monotonically; Figure 10 shows that the proposed algorithms are convergent due to the increase in the number of iterations.

In addition, to further verify the effectiveness of our proposed algorithm, we took the first three training samples from each class on the PolyU palmprint database; Table 4 shows the CPU time of each method.

4.5. Observations and Discussion

The four contributions of the SILR-2DLDGE algorithm are as follows:

• (1). In Table 2 and Table 3, the maximum average recognition rate of the SILR-2DLDGE algorithm on different databases is the highest, which fully shows that our algorithm is optimal and robust. The key reason is that SILR-2DLDGE not only learns a base matrix with low-rank property and local discriminant ability, but also has the advantages of SILR-NMF, which can weaken the disturbance of noise.

• (2). The SILR-2DLDGE algorithm has the advantages of low-rank learning, sparse learning, and incoherent structure learning, similarly to the SILR-NMF algorithm, as well as the advantages of graph embedding, which shows that more discrimination information can be obtained by sparse learning combined with the L1 norm and L * norm. We can see from the curve changes in Figure 3, Figure 4, Figure 5, Figure 7, Figure 8 and Figure 9 that the average recognition rates of this algorithm in three noisy databases and three occlusion databases are higher than that of other algorithms, which fully show that our algorithm has stronger robustness than other algorithms.

• (3). The SILR-2DLDGE algorithm takes less time than the others (Table 4) and learns the sparse transformation matrix to encode the geometric structure of the data and effectively improve the classification accuracy, which can obtain clean data in cases of noise disturbance. At the same time, structural incoherence can ensure that it is easier to separate data points of different classes.

• (4). The maximum average recognition rate of the SILR-2DLDGE algorithm varies among different databases. For example, the SILR-2DLDGE algorithm has the highest average maximum recognition rate without occlusion on the ORL database, whereas it has the lowest average maximum recognition rate with an occlusion of 10 × 10 on the COIL20 database. At the same time, the SILR-2DLDGE algorithm has the highest average maximum recognition rate without noise on the YALE database, whereas it has the lowest average maximum recognition rate with added “salt and pepper” noise with a density of 0.1 on the AR database. The above reasons are due to the small number of ORL and YALE samples, the large AR and PolyU palmprint databases, the presence of glasses and other occlusions in AR, and the presence of illumination and other effects in the PolyU palmprint.

5. Conclusions

This study mainly combined the 2DLPP algorithm with low-rank representation learning and proposed a structurally incoherent low-rank two-dimensional local discriminant graph embedding (SILR-2DLDGE) algorithm based on subspace learning, graph embedding, low-rank sparsity, and structural incoherence. Identification information, local geometry information, low-rank representation information, and structural incoherence existed simultaneously in the proposed algorithm. In addition, the ultimate goal of the proposed algorithm is to make the data points as independent as possible from different classes. In particular, it used low-rank learning L1 norms as constraints to reduce the influence of noise and corruption. The noise and occlusion experiments on six public databases further proved that the proposed SILR-2DLDGE algorithm is more robust than other algorithms. The algorithm is sensitive to parameters. In the future, we will study various parameters and further improve the robustness of the algorithm.

Author Contributions

Conceptualization, H.G. and A.S.; methodology, H.G. and A.S.; software, H.G. and A.S.; validation, H.G. and A.S.; formal analysis, H.G.; investigation, H.G.; resources, H.G.; data curation, H.G.; writing—original draft preparation, H.G. and A.S.; writing—review and editing, H.G. and A.S.; visualization, H.G.; supervision, A.S.; project administration, H.G. and A.S.; funding acquisition, H.G. and A.S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The flow chart of the SILR-2DLDGE method.

Figure 2. Some occlusion images from the three different databases. (a) FERET, (b) ORL, and (c) COIL 100.

Figure 3. The average recognition rates (%) of the eight algorithms in the FERET database vary with the dimension.

Figure 4. The average recognition rates (%) of the eight algorithms in the ORL database varies with the dimension.

Figure 5. The average recognition rates (%) of the eight algorithms in COIL20 database vary with the dimension.

Figure 6. Some the original images and corrosion images on the three different databases. (a) Yale, (b) AR, and (c) PolyU.

Figure 7. The average recognition rates (%) of the eight algorithms in the Yale database vary with the dimension.

Figure 8. The average recognition rates (%) of the eight algorithms in the AR database vary with the dimension.

Figure 9. The average recognition rates (%) of the eight algorithms in the PolyU palmprint database vary with the dimension.

Figure 10. Convergence of the proposed algorithm. (a) FERET databases; (b) AR databases.

Table 1

Parameter values of the SILR-2DLDGE algorithm on different datasets.

Dataset	Parameters of SILR-2DLDGE
	α	β	γ	η
FERET	0.001	0.1	0.3	0.01
ORL	0.01	0.1	0.2	0.001
COIL 100	0.1	10	0.5	100
Yale	0.001	0.01	0.3	0.01
AR	0.001	0.1	0.8	0.01
PolyU	0.01	0.1	0.6	10

Table 2

The best average recognition rates (%) of the eight algorithms in FERET, ORL, and COIL100 databases and the corresponding dimensions (D).

Algorithm	Database
	Without Occlusion			Occlusion Size = 10 × 10
	FERET	ORL	COIL20	FERET	ORL	COIL20
LPP	87.16	89.50	86.69	84.35	87.36	85.90
	25(D)	18(D)	20(D)	26(D)	20(D)	24(D)
2DLPP	89.28	93.11	88.35	86.76	90.55	86.86
	40 × 16(D)	50 × 14(D)	32 × 20(D)	40 × 18(D)	50 × 16(D)	32 × 22(D)
LPP-L1	84.25	86.33	83.45	82.18	85.87	81.65
	28(D)	22(D)	26(D)	30(D)	24(D)	28(D)
2DLPP-L1	85.42	88.36	84.35	83.82	87.58	82.36
	40 × 12(D)	50 × 14(D)	32 × 18(D)	40 × 16(D)	50 ×16(D)	32 × 20(D)
LRR	82.65	86.85	81.20	80.63	85.25	79.95
	32(D)	16(D)	28(D)	36(D)	18(D)	32(D)
LRMD-SI	88.70	91.35	86.80	86.78	89.80	84.95
	16(D)	18(D)	16(D)	18(D)	18(D)	20(D)
SILR-NMF	92.52	93.95	90.65	90.23	91.68	88.28
	418(D)	26(D)	20(D)	22(D)	30(D)	20(D)
SILR-2DLDGE	94.75	96.62	92.54	91.61	94.86	90.58
	40 × 18(D)	50 × 16(D)	32 × 20(D)	40 × 16(D)	50 × 20(D)	32 × 18(D)

D: dimension.

Table 3

The best average recognition rates (%) of different algorithms and corresponding dimensions (D) on the Yale, AR, and PolyU databases and the corresponding dimensions.

Algorithm	Database
	Without Noise			Noise Density = 0.1
	Yale	AR	PolyU	Yale	AR	PolyU
LPP	90.31	88.93	89.08	89.65	85.43	87.75
	32(D)	40(D)	36(D)	36(D)	40(D)	40(D)
2DLPP	91.18	90.35	91.88	90.56	87.39	88.86
	56 × 26(D)	50 × 24(D)	64 × 30(D)	56 × 30(D)	50 × 28(D)	64 × 36(D)
LPP-L1	91.87	89.36	90.56	89.08	86.45	88.82
	36(D)	45(D)	38(D)	32(D)	45(D)	36(D)
2DLPP-L1	92.98	90.16	91.08	90.66	88.45	89.56
	56 × 20(D)	50 × 18(D)	64 × 24(D)	56 × 26(D)	50 × 22(D)	64 × 28(D)
LRR	87.65	82.96	83.25	85.36	80.40	82.82
	36(D)	18(D)	16(D)	40(D)	20(D)	20(D)
LRMD-SI	91.80	86.45	88.23	89.08	83.82	85.66
	20(D)	20(D)	18(D)	22(D)	28(D)	26(D)
SILR-NMF	92.51	90.69	91.48	90.65	87.02	89.54
	12(D)	16(D)	10(D)	14(D)	14(D)	12(D)
SILR-2DLDGE	95.25	93.60	94.78	93.65	90.28	92.15
	56 × 12(D)	50 × 14(D)	64 × 10(D)	56 × 16(D)	50 × 12(D)	64 × 16(D)

D: dimension.

Table 4

The average CPU time consumed of the eight algorithms in the PolyU palmprint database.

PolyUALGORITHM		LPP	2DLPP	$$\mathbf{LPP}- {\mathit{L}}_1$$	$$2\mathbf{DLPP}- {\mathit{L}}_1$$	LRR	LRMD-SI	SILR-NMF	SILR-2DLDGE
Noise density = 0.1	(S)	0.981	0.818	1.236	0.852	0.765	0.783	0.706	0.652

References

1. Qian, J.; Zhu, S.; Wong, W.K.; Zhang, H.; Lai, Z.; Yang, J. Dual robust regression for pattern classifification. Inf. Sci.; 2021; 546, pp. 1014-1029. [DOI: https://dx.doi.org/10.1016/j.ins.2020.09.062]

2. Niu, B.; Yang, Q.; Shiu, S.C.K.; Pal, S.K. Two-dimensional Laplacianfaces algorithm for face recognition. Pattern Recognit.; 2008; 41, pp. 3237-3243. [DOI: https://dx.doi.org/10.1016/j.patcog.2007.12.001]

3. Hu, D.; Feng, G.; Zhou, Z. Two-dimensional locality preserving projections (2DLPP) with its application to palmprint recognition. Pattern Recognit.; 2007; 40, pp. 339-342. [DOI: https://dx.doi.org/10.1016/j.patcog.2006.06.022]

4. Wan, M.; Chen, X.; Zhan, T.; Xu, C.; Yang, G.; Zhou, H. Sparse Fuzzy Two-Dimensional Discriminant Local Preserving Projection (SF2DDLPP) for Robust Image Feature Extraction. Inf. Sci.; 2021; 563, pp. 1-15. [DOI: https://dx.doi.org/10.1016/j.ins.2021.02.006]

5. Zheng, H.; Wang, R.; Ji, W.; Zong, M.; Wong, W.K.; Lai, Z.; Lv, H. Discriminative deep multi-task learning for facial expression recognition. Inf. Sci.; 2020; 533, pp. 60-71. [DOI: https://dx.doi.org/10.1016/j.ins.2020.04.041]

6. Shen, F.; Zhou, X.; Yu, J.; Yang, Y.; Liu, L.; Shen, H.T. Scalable Zero-Shot Learning via Binary Visual-Semantic Embeddings. IEEE Trans. Image Process.; 2019; 28, pp. 3662-3674. [DOI: https://dx.doi.org/10.1109/TIP.2019.2899987] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30794175]

7. Zhao, C.; Wang, X.; Zuo, W.; Shen, F.; Shao, L.; Miao, D. Similarity learning with joint transfer constraints for person re-identification. Pattern Recognit.; 2020; 97, 107014. [DOI: https://dx.doi.org/10.1016/j.patcog.2019.107014]

8. Zhou, J.; Lai, Z.; Miao, D.; Gao, C.; Yue, X. Multigranulation rough-fuzzy clustering based on shadowed sets. Inf. Sci.; 2020; 507, pp. 553-573. [DOI: https://dx.doi.org/10.1016/j.ins.2018.05.053]

9. Zhang, Z.; Zhang, Y.; Liu, G.; Tang, J.; Yan, S.; Wang, M. Joint Label Prediction Based Semi-Supervised Adaptive Concept Factorization for Robust Data Representation. IEEE Trans. Knowl. Data Eng.; 2019; 32, pp. 952-970. [DOI: https://dx.doi.org/10.1109/TKDE.2019.2893956]

10. Liu, Z.; Lai, Z.; Ou, W.; Zhang, K.; Zheng, R. Structured optimal graph based sparse feature extraction for semi-supervised learning. Signal Process; 2020; 170, 107456. [DOI: https://dx.doi.org/10.1016/j.sigpro.2020.107456]

11. Zhi, R.; Ruan, Q. Facial expression recognition based on two-dimensional discriminant locality preserving projections. Neurocomputing; 2008; 71, pp. 1730-1734. [DOI: https://dx.doi.org/10.1016/j.neucom.2007.12.002]

12. Wan, M.; Lai, Z.; Shao, J.; Jin, Z. Two-dimensional local graph embedding discriminant analysis (2DLGEDA) with its application to face and Palm Biometrics. Neurocomputing; 2009; 73, pp. 197-203. [DOI: https://dx.doi.org/10.1016/j.neucom.2009.07.015]

13. Jove, E.; Casteleiro-Roca, J.-L.; Quintián, H.; Méndez-Pérez, J.-A.; Calvo-Rolle, J.-L. A new method for anomaly detection based on non-convex boundaries with random two-dimensional projections. Inf. Fusion; 2021; 65, pp. 50-57. [DOI: https://dx.doi.org/10.1016/j.inffus.2020.08.011]

14. Meng, D.; Zhao, Q.; Xu, Z. Improve robustness of sparse PCA by L1-norm maximization. Pattern Recognit.; 2012; 45, pp. 487-497. [DOI: https://dx.doi.org/10.1016/j.patcog.2011.07.009]

15. Kwak, N. Principal Component Analysis Based on L1-Norm Maximization. IEEE Trans. Pattern Anal. Mach. Intell.; 2008; 30, pp. 1672-1680. [DOI: https://dx.doi.org/10.1109/TPAMI.2008.114]

16. Ding, C.; Zhou, D.; He, X.; Zha, H. R1-PCA: Rotational invariant L1-norm principal component analysis for robust subspace factorization. Proceedings of the 23rd International Conference on Machine Learning; Pittsburgh, PA, USA, 25–29 June 2006; pp. 281-288.

17. Pang, Y.; Yuan, Y. Outlier-resisting graph embedding. Neurocomputing; 2010; 73, pp. 968-974. [DOI: https://dx.doi.org/10.1016/j.neucom.2009.08.020]

18. Zhao, H.X.; Xing, H.J.; Wang, X.Z.; Chen, J.F. L1-norm-based 2DLPP. Proceedings of the 2011 Chinese Control and Decision Conference (CCDC); Mianyang, China, 23–25 May 2011; pp. 1259-1264.

19. Liu, G.; Lin, Z.; Yu, Y. Robust subspace segmentation by low-rank representation. Proceedings of the 27th International Conference on International Conference on Machine Learning; Haifa, Israel, 21–24 June 2010; pp. 663-670.

20. Liu, G.; Lin, Z.; Yan, S.; Sun, J.; Yu, Y.; Ma, Y. Robust Recovery of Subspace Structures by Low-Rank Representation. IEEE Trans. Pattern Anal. Mach. Intell.; 2012; 35, pp. 171-184. [DOI: https://dx.doi.org/10.1109/TPAMI.2012.88]

21. Wright, J.; Ganesh, A.; Rao, S.; Peng, Y.; Ma, Y. Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization. Proceedings of the 22nd International Conference on Neural Information Processing Systems; Whistler, BC, Canada, 7–10 December 2009; pp. 1-9.

22. Candès, E.J.; Li, X.; Ma, Y.; Wright, J. Robust principal component analysis. J. ACM; 2009; 58, pp. 1-17. [DOI: https://dx.doi.org/10.1145/1970392.1970395]

23. Liu, J.; Chen, Y.; Zhang, J.; Xu, Z. Enhancing Low-Rank Subspace Clustering by Manifold Regularization. IEEE Trans. Image Process.; 2014; 23, pp. 4022-4030. [DOI: https://dx.doi.org/10.1109/TIP.2014.2343458]

24. Zhuang, L.; Gao, H.; Lin, Z.; Ma, Y.; Zhang, X.; Yu, N. Non-negative low-rank and sparse graph for semi-supervised learning. Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition; Providence, RI, USA, 16–21 June 2012; pp. 2328-2335.

25. Lu, Y.; Yuan, C.; Zhu, W.; Li, X. Structurally Incoherent Low-Rank Nonnegative Matrix Factorization for Image Classification. IEEE Trans. Image Process.; 2018; 27, pp. 5248-5260. [DOI: https://dx.doi.org/10.1109/TIP.2018.2855433]

26. Wei, C.-P.; Chen, C.-F.; Wang, Y.-C.F. Robust face recognition with structurally incoherent low-rank matrix decomposition. IEEE Trans. Image Process.; 2014; 23, pp. 3294-3307. [DOI: https://dx.doi.org/10.1109/TIP.2014.2329451]

27. Wright, J.; Yang, A.Y.; Ganesh, A.; Sastry, S.S.; Ma, Y. Robust face recognition via sparse representation. IEEE Trans. Pattern Anal. Mach. Intell.; 2009; 31, pp. 210-227. [DOI: https://dx.doi.org/10.1109/TPAMI.2008.79]

28. Lin, Z.; Chen, M.; Ma, Y. The augmented Lagrange multiplier algorithm for exact recovery of corrupted low-rank matrices. arXiv; 2009; arXiv: 1009.5055

29. Boyd, S.; Parikh, N.; Chu, E.; Peleato, B.; Eckstein, J. Distributed optimization and statistical learning via the alternating direction algorithm of multipliers. Found. Trends Mach. Learn.; 2011; 3, pp. 1-122. [DOI: https://dx.doi.org/10.1561/2200000016]

30. Mazumder, R.; Hastie, T.; Tibshirani, R. Spectral Regularization Algorithms for Learning Large Incomplete Matrices. J. Mach. Learn. Res.; 2010; 11, pp. 2287-2322.

31. Drineas, P.; Frieze, A.; Kannan, R.; Vempala, S.; Vinay, V. Clustering Large Graphs via the Singular Value Decomposition. Mach. Learn.; 2004; 56, pp. 9-33. [DOI: https://dx.doi.org/10.1023/B:MACH.0000033113.59016.96]

32. Hosseini-Asl, E.; Zurada, J.M.; Nasraoui, O. Deep Learning of Part-Based Representation of Data Using Sparse Autoencoders With Nonnegativity Constraints. IEEE Trans. Neural Netw. Learn. Syst.; 2015; 27, pp. 2486-2498. [DOI: https://dx.doi.org/10.1109/TNNLS.2015.2479223]

33. Boyd, S.P.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004.

34. Kim, E.; Lee, M.; Oh, S. Elastic-net regularization of singular values for robust subspace learning. Proceedings of the 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR); Boston, MA, USA, 7–12 June 2015; pp. 915-923.

35. He, X.; Niyogi, P. Locality Preserving Projections. Proceedings of the Advances in Neural Information Processing System 16; Vancouver, BC, Canada, 8–13 December 2003.

36. He, X.; Yan, S.; Hu, Y.; Niyogi, P.; Zhang, H.J. Face recognition using laplacianfaces. IEEE Trans. Pattern Anal. Mach. Intell.; 2005; 27, pp. 328-340.

37. Phillips, P.J.; Moon, H.; Rizvi, S.A.; Rauss, P.J. The FERET evaluation algorithmology for face-recognition algorithms. IEEE Trans. Pattern Anal. Mach. Intell.; 2000; 22, pp. 1090-1104. [DOI: https://dx.doi.org/10.1109/34.879790]

38. Samaria, F.S.; Harter, A.C. Parameterisation of a stochastic model for human face identifification. Proceedings of the 1994 IEEE Workshop on Applications of Computer Vision; Sarasota, FL, USA, 5–7 December 1994; pp. 138-142.

39. Nene, S.A.; Nayar, S.K.; Murase, H. Columbia Object Image Library (COIL-100); Columbia University: New York, NY, USA, 1996.

40. Yang, J.; Zhang, D.; Frangi, A.F.; Yang, J.Y. Two-dimensional PCA: A new approach to appearance-based face representation and recognition. IEEE Trans. Pattern Anal. Mach. Intell.; 2004; 26, pp. 131-137. [DOI: https://dx.doi.org/10.1109/TPAMI.2004.1261097]

41. Martinez, A.M.; Benavente, R. The AR Face Database; Centrede Visió per Computador, Univ. Autonoma de Barcelona: Bellaterra, Barcelona, 1998.

42. Liu, Z.; Ou, W.; Lu, W.; Wang, L. Discriminative feature extraction based on sparse and low-rank representation. Neurocomputing; 2019; 362, pp. 129-138. [DOI: https://dx.doi.org/10.1016/j.neucom.2019.06.073]

43. Yang, J.; Zhang, D.; Yang, J.; Niu, B. Globally maximizing, locally minimizing: Unsupervised discriminant projection with applications to face and palm biometrics. IEEE Trans. Pattern Anal. Mach. Intell.; 2007; 29, pp. 650-664. [DOI: https://dx.doi.org/10.1109/TPAMI.2007.1008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17299222]