Relaxed 2-D Principal Component Analysis by $L_p$ Norm for Face Recognition

TL;DR

This paper introduces a relaxed 2D PCA method that incorporates label information and $L_p$ norms to improve face recognition accuracy, outperforming existing techniques.

Contribution

It proposes a novel R2DPCA approach using label-based relaxation vectors and $L_p$ norms, enhancing recognition performance over traditional 2DPCA methods.

Findings

R2DPCA achieves higher recognition rates than state-of-the-art methods.

The method demonstrates high generalization ability on practical face datasets.

Optimal $L_p$-norms are effectively selected within a reasonable range.

Abstract

A relaxed two dimensional principal component analysis (R2DPCA) approach is proposed for face recognition. Different to the 2DPCA, 2DPCA-$L_1$ and G2DPCA, the R2DPCA utilizes the label information (if known) of training samples to calculate a relaxation vector and presents a weight to each subset of training data. A new relaxed scatter matrix is defined and the computed projection axes are able to increase the accuracy of face recognition. The optimal $L_p$-norms are selected in a reasonable range. Numerical experiments on practical face databased indicate that the R2DPCA has high generalization ability and can achieve a higher recognition rate than state-of-the-art methods.

Tables

Table 1. Table 1 : Procedures

Algorithm	2DPCA	2DPCA-$L_1$
Procedure	$[\text{W},\text{D}]=\text{eig}\left(\frac{1}{n}\sum _{i=1}^n{\text{X}}_i^T{\text{X}}_i\right)$	${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T\text{sign}\left({\text{X}}_i{\text{w}}^k\right)$
		${\text{w}}^{k+1}=\frac{{\text{v}}^k}{{\Vert {\text{v}}^k\Vert }_2}$
Algorithm	2DPCA$L_1$-S	G2DPCA
Procedure	${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T\text{sign}\left({\text{X}}_i{\text{w}}^k\right)$	${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T\left[{\left|{\text{X}}_i{\text{w}}^k\right|}^{s-1}\circ \text{sign}\left({\text{X}}_i{\text{w}}^k\right)\right]$
	$u_i^k=v_i^k\frac{\left|w_i^k\right|}{\lambda +\left|w_i^k\right|},i=1,2,\dots ,w$	if $p\ge 1$    run (9)
	${\text{w}}^{k+1}=\frac{{\text{u}}^k}{{\Vert {\text{u}}^k\Vert }_2}$	if    run (12)

Table 2. Table 2 : Classification accuracies of five algorithms on faces95

$Algorithms$	$Optimalparameters$	$Accuracy$
2DPCA	$-$	$0.8729$
2DPCA-$L_1$	$-$	$0.8708$
2DPCA$L_1$-S	$\rho =-0.5$	$0.8785$
G2DPCA	$s=2.7,p=2.2$	$0.9451$
R2DPCA $\left(\gamma =0\right)$	$s=1,p=2.2$	0.9493

Table 3. Table 3 : Classification accuracies of five algorithms on Color Feret

$Algorithms$	$Optimalparameters$	$Accuracy$
2DPCA	$-$	$0.5982$
2DPCA-$L_1$	$-$	$0.5985$
2DPCA$L_1$-S	$\rho =-0.3$	$0.6236$
G2DPCA	$s=2.8,p=2.6$	$0.6918$
R2DPCA $\left(\gamma =0\right)$	$s=3,p=2.2$	0.7085

Table 4. Table 4 : Classification accuracies according to different γ 𝛾 \gamma on faces95.

Database	$\gamma$	Optimal parameters	Accuracy
faces95	$0:1/4:1$	$s=1.1,p=2.2$	$0.8861$
color FERET	$0:1/4:3/4$	$s=3,p=2.2$	$0.7673$
	$1$	$s=3,p=2.2$	$0.7733$

Table 5. Table 5 : Classification accuracies according to different γ 𝛾 \gamma on Gray Feret.

$\gamma$	Optimal parameters	Accuracy
$0$	$s=1.8,p=1.7$	$0.6075$
$1/4$	$s=1.9,p=1.7$	$0.6075$
$1/2$	$s=2.3,p=1.6$	$0.6112$
$3/4$	$s=2.4,p=1.6$	$0.6088$
$1$	$s=2.6,p=1.8$	$0.5837$

Table 6. Table 6 : Restarted alternating searching path on Face95.

$initialpoint$	$path$	$searchingAccuracy$
$s=1,p=2.2$	$(1,2.2)$	$0.9347$
$s=1.3,p=0.9$	$(1.3,0.9)\to (1.3,2.1)\to (1.8,2.1)$	0.9375
$s=1.6,p=2.8$	$(1.6,2.8)\to (1.3,2.8)$	$0.9306$
$s=2.7,p=1.2$	$(2.7,1.2)\to (2.7,2.2)\to (1,2.2)\to (1,2.1)\to (1.8,2.1)$	0.9375

Table 7. Table 7 : Restarted alternating searching path on Color Feret.

$initialpoint$	$path$	$searchingAccuracy$
$s=1,p=2.4$	$(1,2.4)\to (2.7,2.4)\to (2.7,2.2)$	0.7158
$s=1.7,p=0.9$	$(1.7,0.9)\to (1.7,2.3)\to (1.6,2.3)\to (1.6,2.6)\to (2.5,2.6)\to (2.5,2.1)$	$0.7085$
$s=1.6,p=1.7$	$(1.6,1.7)\to (2.5,2.7)\to (2.5,2.2)\to (2.5,2.2)\to (2.7,2.2)$	0.7158
$s=1.2,p=1.2$	$(1.2,1.2)\to (1.2,2.3)$	$0.6939$

Table 8. Table 8 : classification accuracies of six algorithms on Gray FERET

$Algorithms$	$Optimalparameters$	$Accuracy$
2DPCA	$-$	$0.4225$
R2DPCA $\left(s=p=2\right)$	$-$	0.6387
2DPCA-$L_1$	$-$	$0.4225$
R2DPCA$\left(s=1,p=2\right)$	$-$	0.6123
2DPCA$L_1$-S	$\rho =-0.6$	$0.4938$
R2DPCA $\left(s=1\right)$	$\rho =-2$	0.5988

Equations

$${\text{X}}_i^{deflated}={\text{X}}_i(\text{I}-\text{W}{\text{W}}^T),i=1,2,\dots ,n.$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }\sum _{i=1}^n\Vert {\text{X}}_i\text{w}{\Vert }_2^2,s.t.\Vert \text{w}{\Vert }_2^2=1.$$

$${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T[|{\text{X}}_i{\text{w}}^k|\circ \text{sign}({\text{X}}_i{\text{w}}^k)],$$

$${\text{u}}^k=|{\text{v}}^k|\circ \text{sign}({\text{v}}^k),$$

$${\text{w}}^{k+1}=\frac{{\text{u}}^k}{\Vert {\text{u}}^k{\Vert }_2},$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }\sum _{i=1}^n\Vert {\text{X}}_i\text{w}{\Vert }_1,s.t.\Vert \text{w}{\Vert }_2^2=1.$$

$${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T\text{sign}({\text{X}}_i{\text{w}}^k),$$

$${\text{w}}^{k+1}=\frac{{\text{v}}^k}{\Vert {\text{v}}^k{\Vert }_2},$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }\sum _{i=1}^n\Vert {\text{X}}_i\text{w}{\Vert }_1,s.t.\Vert \text{w}{\Vert }_1\le c,\Vert \text{w}{\Vert }_2^2=1,$$

$${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T\text{sign}({\text{X}}_i{\text{w}}^k),$$

$$u_i^k=v_i^k\frac{|w_i^k|}{\lambda +|w_i^k|},i=1,2,\dots ,w,$$

$${\text{w}}^{k+1}=\frac{{\text{u}}^k}{\Vert {\text{u}}^k{\Vert }_2},$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }\sum _{i=1}^n\Vert {\text{X}}_i\text{w}{\Vert }_s^s,s.t.\Vert \text{w}{\Vert }_p^p=1$$

$${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T[|{\text{X}}_i{\text{w}}^k{|}^{s-1}\circ \text{sign}({\text{X}}_i{\text{w}}^k)],$$

$${\text{u}}^k=|{\text{v}}^k{|}^{q-1}\circ \text{sign}({\text{v}}^k),$$

$${\text{w}}^{k+1}=\frac{{\text{u}}^k}{\Vert {\text{u}}^k{\Vert }_p}.$$

$$j=\arg {\max }_{i\in [1,w]}|v_i^k|,$$

$$w_i^{k+1}=\{\begin{array}{r}\text{sign}(v_j^k),i=j,\\ 0,i\ne j,\end{array}$$

$${\text{w}}^{k+1}=\text{sign}({\text{v}}^k).$$

$${\text{v}}^k=\sum _{i=1}^n{\text{X}}_i^T[|{\text{X}}_i{\text{w}}^k{|}^{s-1}\circ \text{sign}({\text{X}}_i{\text{w}}^k)],$$

$${\text{u}}^k=|{\text{w}}^k{|}^{2-p}\circ {\text{v}}^k,$$

$${\text{w}}^{k+1}=\frac{{\text{u}}^k}{\Vert {\text{u}}^k{\Vert }_p}.$$

$${\text{G}}_t=\frac{1}{n}\sum _{i=1}^n{\text{X}}_i^T{\text{X}}_i,$$

$$[\text{W},\text{D}]=\text{eig}({\text{G}}_t),$$

$$\{{\text{w}}_1,{\text{w}}_2,\cdots ,{\text{w}}_r\}=\arg \max \sum _{i=1}^n{\text{w}}_i^T{\text{G}}_t{\text{w}}_i,$$

$$\textbf{w}_{i}^{T}\textbf{w}_{j}=\left\{\begin{array}[]{l}1,~{}i=j,\\ 0,~{}i\neq j,\end{array}~{}s,t=1,\cdots,r.\right.$$

$${\text{X}}_1^1,\cdots ,{\text{X}}_{n_1}^1|{\text{X}}_1^2,\cdots ,{\text{X}}_{n_2}^2|\cdots |{\text{X}}_1^m,\cdots ,{\text{X}}_{n_m}^m,$$

$${\text{M}}_j=\frac{1}{n_j}\sum _{i=1}^{n_j}{\text{X}}_i^j\in {\mathbb{R}}^{h\times w},$$

$${\text{C}}_j=\frac{1}{n_j}\sum _{i=1}^{n_j}({\text{X}}_i^j-{\text{M}}_j{)}^T({\text{X}}_i^j-{\text{M}}_j)\in {\mathbb{R}}^{w\times w},$$

$$\text{v}=[v_1,\cdots ,v_m{]}^T\in {\mathbb{R}}^m,$$

$$v_j=\frac{f({\lambda }_{\max }({\text{C}}_j))}{{\sum }_{i=1}^{m}f({\lambda }_{\max }({\text{C}}_i))},$$

$$\text{M}=\frac{1}{n}\sum _{i=1}^n{\text{X}}_i=\frac{1}{n}\sum _{j=1}^m\sum _{i=1}^{n_j}{\text{X}}_i^j.$$

$$J(\text{w})=\gamma \text{G}+(1-\gamma )\tilde{\text{G}},$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }J(\text{w}),s.t.\Vert \text{w}{\Vert }_p^p=1,$$

$$\underset{\text{w}\in {\mathbb{R}}^w}{\max }\sum _{j=1}^m\sum _{i=1}^{n_j}\Vert \frac{v_j}{n_j}({\text{X}}_i^j-\text{M})\text{w}{\Vert }_s^s,s.t.\Vert \text{w}{\Vert }_p^p=1.$$

$$f_j=\gamma \sum _{i=1}^n\Vert ({\text{X}}_i-\text{M}{)}^{deflated}{\text{w}}_j{\Vert }_s^s+(1-\gamma )\sum _{j=1}^m\sum _{i=1}^{n_j}\Vert \frac{v_j}{n_j}({\text{X}}_i^j-\text{M}{)}^{deflated}{\text{w}}_j{\Vert }_s^s.$$

Full text

Key words. Face recognition; G2DPCA; Relaxed 2DPCA; Optimal algorithms; Alternating direction method

1 Introduction

The principal component analysis (PCA) [1, 2], has become one of the most powerful approaches of face recognition [3, 4, 5, 6, 7]. Recently, many robust PCA (RPCA) algorithms are proposed with improving the quadratic formulation, which renders PCA vulnerable to noises, into $L_{1}$-norm on the objection function, e.g., $L_{1}$-PCA [8], $R_{1}$-PCA [9], and PCA-$L_{1}$ [10]. Meanwhile, sparsity is also introduced into PCA algorithms, resulting in a series of sparse PCA (SPCA) algorithms [11, 12, 13, 14]. A newly proposed robust SPCA (RSPCA) [15] further applies $L_{1}$-norm both in objective and constraint functions of PCA, inheriting the merits of robustness and sparsity. Observing that $L_{2}$-, $L_{1}$-, and $L_{0}$-norms are all special $L_{p}$-norm, it is natural to impose $L_{p}$-norm on the objection or/and constraint functions, straightforwardly; see PCA-$L_{p}$ [16] and generalized PCA (GPCA) [17] for instance.

To preserve the spatial structure of face images, two dimensional PCA (2DPCA), proposed by Yang et al. [18], represents face images with two dimensional matrices rather than one dimensional vectors. The computational problems bases on 2DPCA are of much smaller scale than those based on traditional PCA, and the difficulties caused by low rank are also avoided. This image-as-matrix method offers insights for improving above RSPCA, PCA-$L_{p}$, GPCA, etc. As typical examples, the $L_{1}$-norm-based 2DPCA (2DPCA-$L_{1}$) [19] and 2DPCA-$L_{1}$ with sparsity (2DPCA$L_{1}$-S) [20] are improvements of PCA-$L_{1}$ and RSPCA, respectively, and the generalized 2DPCA (G2DPCA) [21] imposes $L_{p}$-norm on both objective and constraint functions of 2DPCA. Recently, the quaternion 2DPCA is proposed in [22] and applied to color face recognition, where the red, green and blue channels of a color image is encoded as three imaginary parts of a pure quaternion matrix. To arm the quaternion 2DPCA with the generalization ability, Zhao, Jia and Gong [23] proposed the sample-relaxed quaternion 2DPCA with applying the label information (if known) of training samples. The structure-preserving algorithms of quaternion eigenvalue decomposition and singular value decomposition can be found in [24, 25, 26, 27, 28, 29, 30, 31]. More applications of the quaternion representation and structure-preserving methods to color image processing can be found in [32] and [33].

Both PCA and 2DPCA are unsupervised methods and omit the potential or known label information of samples. They are often applied to the training set and thus the computed projections will maximize the scatter of projected training samples. That means the scatter of projected testing samples are not surely optimal, and certainly, so are the whole (training and testing) projected samples. Inspired by this observation, we proposed a new relaxation two-dimensional principal component analysis (R2DPCA) in this paper. R2DPCA sufficiently utilizes the labels (if known) of training samples, and can enhance the total scatter of whole projected samples. This approach is a generalization of G2DPCA [21], and will reduce to G2DPCA if the label information is unknown or unused.

The rest of this paper is organized as follows. In Section 2, we recall robust and sparse 2DPCA algorithms. In Section 3, we present a new relaxed two dimensional principal component analysis (R2DPCA) approach for face recognition. In Section 4, we compare the R2DPCA with the state-to-the-art approaches, and indicate the efficiencies of the R2DPCA . In Section 5, we sum up the contribution of this paper.

2 Robust and sparse 2DPCA algorithms

In this section, we recall 2DPCA, 2DPCA-$L_{1}$, 2DPCA$L_{1}$-S, and G2DPCA algorithms in the form of computing the first projection vector. In fact, after obtaining first $j$ projection vectors $\textbf{W}=[\textbf{w}_{1},\textbf{w}_{2},\ldots,\textbf{w}_{j}]$, the $(j+1)$-th projection vector $\textbf{w}_{j+1}$ can be calculated similarly on deflated samples [34]:

[TABLE]

2.1 2DPCA

Suppose that there are $n$ training images samples $\textbf{X}_{1},\textbf{X}_{2},\ldots,\textbf{X}_{n}\in\mathbb{R}^{h\times w}$, where $h$ and $w$ denote the height and width of images, respectively. We assume that these samples are mean-centered, i.e., $\frac{1}{n}\sum^{n}\nolimits_{i=1}\textbf{X}_{i}=0$; otherwise, we will replace $\textbf{X}_{i}$ by $\textbf{X}_{i}-\frac{1}{n}\sum^{n}\nolimits_{i=1}\textbf{X}_{i}$.

2DPCA [18] finds its first projection vector $\textbf{w}\in\mathbb{R}^{w}$ by solving the optimization problem with equality constraints:

[TABLE]

The projection vector w could be calculated by the iterative algorithm:

[TABLE]

where sign$(\cdot)$ denotes the sign function. The projection vector w can also be obtained by calculating the eigen decomposition of a covariance matrix and selecting the eigenvector corresponding to the largest eigenvalue. See Remark 2.2 for more details.

2.2 2DPCA-$L_{1}$

2DPCA-$L_{1}$ [19] finds its first projection vector $\textbf{w}\in\mathbb{R}^{w}$ by solving the optimization problem with equality constraints:

[TABLE]

The projection vector w could be calculated by the iterative algorithm:

[TABLE]

where $\textbf{w}^{k}$ is the projection vector at the $k$-th step. Notice that 2DPCA-$L_{1}$ could be formulated by replacing the $L_{2}$-norm in objective function of 2DPCA with $L_{1}$-norm.

2.3 2DPCA$L_{1}$-S

2DPCA$L_{1}$-S [20] finds its first projection vector $\textbf{w}\in\mathbb{R}^{w}$ by solving the optimization problem with equality and inequality constraints:

[TABLE]

where $c$ is a positive constant. The projection vector w could be calculated by the iterative algorithm:

[TABLE]

where $u_{i}^{k}$, $v_{i}^{k}$, and $w_{i}^{k}$ are the $i$th elements of vectors $\textbf{u}^{k},\textbf{v}^{k},$ and $\textbf{w}^{k}\in\mathbb{R}^{w}$, respectively. In equation (7b), $\lambda$ is a positive scalar which serves as a tuning parameter. When $\lambda$ is set to be zero, 2DPCA$L_{1}$-S reduces to 2DPCA-$L_{1}$. Notice that 2DPCA$L_{1}$-S could be formulated by imposing $L_{1}$-norm on objective and constraint functions of 2DPCA.

2.4 G2DPCA

G2DPCA [21] finds its first projection vector $\textbf{w}\in\mathbb{R}^{w}$ by solving the optimization problem with equality constraints:

[TABLE]

where $s\geq 1$ and $p>0$. The projection vector w can be updated in two different ways, depending on the value $p$.

• Case 1: If $p\geq 1$,

[TABLE]

where $q$ satisfies $1/p+1/q=1$, $\circ$ denotes the Hadamard product, i.e., the element-wise product between two vectors. Especially, if $p=1$, $\textbf{w}^{k+1}$ can be computed by

[TABLE]

wherein $v_{i}^{k}$ is the $i$-th value of $\textbf{v}^{k}$; if $p=\infty$, $\textbf{w}^{k+1}$ can be computed by

[TABLE]

• Case 2: If $0<p<1$,

[TABLE]

Notice that G2DPCA could be formulated by generalizing $L_{2}$-norm in objective and constraint functions of 2DPCA to $L_{s}$-norm and $L_{p}$-norm, respectively.

Remark 2.1**.**

When $s=p=2$, the projection method (9) reduces to (3).

Remark 2.2**.**

By the eigenvalue decomposition method, 2DPCA can select a set of projection axes $\{\textbf{w}_{1},\textbf{w}_{2},\cdots,\textbf{w}_{r}\}$ in one step, without selecting only one optimal projection axis each step. These projection axes are chosen as eigenvectors of a covariant matrix corresponding to first $r$ largest eigenvalues:

[TABLE]

where $\textbf{G}_{t}$ is the covariance matrix of training samples, $\textbf{W}=[\textbf{w}_{1},\textbf{w}_{2},\cdots,\textbf{w}_{r}]\in R^{w\times r}$ is a matrix with unitary column vectors, D is a diagonal matrix consists of first $r$ largest eigenvalues. Since $G_{t}$ is symmetric and positive semi-definite, the diagonal elements of $D$ are nonnegative and the projection axes

[TABLE]

are orthogonal to each other, i.e.,

[TABLE]

Remark 2.3**.**

We sum the procedures of above four methods in Table 1. Their relationship is also indicated in Fig 2.

3 The relaxed 2DPCA by $L_{p}$-norm

In this section, we introduce a relaxed two-dimensional principal component analysis (R2DPCA) method by $L_{p}$-norm. R2DPCA includes three parts: relaxation vector generation, objective function relaxation, and projection relaxation.

3.1 Relaxation vector

Suppose that training samples $\textbf{X}_{1},\textbf{X}_{2},...,\textbf{X}_{n}\in\mathbb{R}^{h\times w}$ can be partitioned into $m$ classes and each class contains $n_{j}$ samples:

[TABLE]

where $\textbf{X}_{i}^{j}$ denotes the $i$-th sample of the $j$-th class, $i=1,\ldots,n_{j}$, $j=1,\ldots,m$. Define the mean of training samples from the $j$-th class as

[TABLE]

and the $j$-th within-class covariance matrix of the training set as

[TABLE]

where $j=1,\ldots,m$, $\sum_{j=1}^{m}n_{j}=n$ and $i=1,\ldots,n_{j}$.

The within-class covariance matrix $\textbf{C}_{j}$ is a symmetric and positive semi-definite matrix. Its maximal eigenvalue, denoted by $\lambda_{\max}(\textbf{C}_{j})$, represents the variance of training samples $\textbf{X}_{1}^{j},\ldots,\textbf{X}_{n_{j}}^{j}$ in the principal component. In general, the larger $\lambda_{\rm max}(\textbf{C}_{j})$ is, the better scattered of the training samples of $j$-th class are. A very small $\lambda_{\rm max}(\textbf{C}_{j})$ indicates that $\textbf{X}_{1}^{j},\ldots,\textbf{X}_{n_{j}}^{j}$ are not well scattered samples to represent the $j$-th class. Extremely, if $\lambda_{\max}(\textbf{C}_{j})=0$ then all of training samples from the $j$-th class are same, and then the contribution of the $j$-th class to the covariance matrix of training set should be controlled by a small factor. To this aim, we define a relaxation vector of training classes,

[TABLE]

where

[TABLE]

is a relaxation factor of the $j$-th class with a function, $f:~{}\mathbb{R}\rightarrow\mathbb{R}^{+}$. A relaxation factor of each training sample of $j$-th class is defined as $v_{j}/n_{j}$. If each training class has only one sample, i.e., $n_{1}=\cdots=n_{m}=1$, then all within-class covariance matrices are zero matrix and $\lambda_{\max}(\textbf{C}_{1})=\cdots=\lambda_{\max}(\textbf{C}_{m})=0.$ In this case, the relaxation factor of each class is same ($v_{j}=1/m$), and so is the factor of each training sample of $j$-th class.

We sum above steps of computing the relaxation vector of training set in Algorithm 3.1.

Algorithm 3.1**.**

Relaxation vector generation**

*  function $\textbf{v}=\texttt{relaxvec}(\textbf{X}_{1},\textbf{X}_{2},\cdots,\textbf{X}_{n},m,w)$ *

*  for $j=1,2,\cdots,m$ do

*     $\textbf{C}_{j}=\texttt{zeros}(w,w);$ *

*     $\textbf{M}_{j}=\frac{1}{n_{j}}(\textbf{X}_{1}^{j}+\cdots+\textbf{X}_{n_{j}}^{j});$ *

*     for $i=1,2,\cdots,n_{j}$ do

*        $\textbf{C}_{j}=\textbf{C}_{j}+{(\textbf{X}_{i}^{j}-\textbf{M}_{j})}^{\prime}*(\textbf{X}_{i}^{j}-\textbf{M}_{j});$ *

*     end *for

*     $\textbf{C}_{j}=\textbf{C}_{j}/n_{j};$ *

*     Compute relaxation vector v defined as in (16) and (17); *

*  end *for

3.2 Objective function relaxation

Let M denote the mean of training samples, i.e.,

[TABLE]

With computed relaxation vector $\textbf{v}=$ $[v_{1},$ $\cdots,$ $v_{m}]^{T}$ in Section 3.1, we define a relaxed criterion as

[TABLE]

where $\gamma\in[0,1]$ is a relaxation parameter, $\textbf{w}\in\mathbb{R}^{w}$ is a unit vector under $L_{p}$ norm, $\textbf{G}:=\sum\limits^{n}_{i=1}\|(\textbf{X}_{i}-\textbf{M})\textbf{w}\|_{s}^{s}$ and $\widetilde{\textbf{G}}:=\sum\limits^{m}_{j=1}\sum\limits^{n_{j}}_{i=1}\|\frac{v_{j}}{n_{j}}(\textbf{X}_{i}^{j}-\textbf{M})\textbf{w}\|_{s}^{s}$. R2DPCA finds its first projection vector $\textbf{w}\in\mathbb{R}^{w}$ by solving the optimization problem with equality constraints:

[TABLE]

where the criterion $J(\textbf{w})$ is defined as in (18). Notice that (19) reduces to (8) if $\gamma=1$, and thus, the first projection vector of R2DPCA is the same as that of G2DPCA. When $\gamma=0$, (19) is simplified as

[TABLE]

If first $j$ projection vectors $\textbf{W}=[\textbf{w}_{1},\textbf{w}_{2},\ldots,\textbf{w}_{j}]$ have been obtained, the $(j+1)$-th projection vector $\textbf{w}_{j+1}$ can be calculated similarly on the deflated samples, defined as in (1). From each iterative step, we also obtain a maximized objective function value corresponding to $\textbf{w}_{j}$,

[TABLE]

With the relaxed criterion defined in (18), first $r$ optimal projection vectors of R2DPCA solve the optimal problem with equality constraints:

[TABLE]

We propose Algorithm 3.2 to compute first $r$ optimal projection vectors, $\textbf{W}=[\textbf{w}_{1},\cdots,\textbf{w}_{r}]$, and corresponding optimal objective function values, $\textbf{D}={\rm diag}(f_{1},\ldots,f_{r})$.

Algorithm 3.2**.**

R2DPCA**

0:*  $\textbf{X}_{1},\textbf{X}_{2},\cdots,\textbf{X}_{n},s\in[1,\infty),p\in(0,\infty],r,m,w,\gamma\in[0,1],n_{1},\cdots,n_{m},tol.$ *

0:*  $\textbf{W}=[\textbf{w}_{1},\ldots,\textbf{w}_{r}],\textbf{D}=\texttt{diag}(f_{1},\ldots,f_{r}).$ *

*  Initialize $\textbf{W}=[\ ]$, $\textbf{D}=[\ ]$. *

*  v=relaxvec($\textbf{X}_{1},\textbf{X}_{2},\cdots,\textbf{X}_{n},m,w$). *

*  Homogenize training samples. *

*  for $t=1,2,\cdots,r$ do

*     Initialize $k=0$, $\delta=1$, arbitrary $\textbf{w}^{0}$ with $\parallel\textbf{w}^{0}\parallel_{p}=1$. *

*     $f_{0}=\gamma\sum\limits^{n}_{i=1}\|\textbf{X}_{i}\textbf{w}^{0}\|_{s}^{s}+(1-\gamma)\sum\limits^{m}_{j=1}\sum\limits^{n_{j}}_{j=1}\|\frac{v_{j}}{n_{j}}\textbf{X}_{i}^{j}\textbf{w}^{0}\|_{s}^{s}.$ *

*     while $\delta>tol$ do

*        $\textbf{v}^{k}=\gamma\sum\limits^{n}_{i=1}\textbf{X}_{i}^{T}[|\textbf{X}_{i}\textbf{w}^{k}|^{s-1}\circ\texttt{sign}(\textbf{X}_{i}\textbf{w}^{k})]+(1-\gamma)\sum\limits^{m}_{j=1}\sum\limits^{n_{j}}_{j=1}(\frac{v_{j}}{n_{j}}\textbf{X}_{i}^{j})^{T}[|\frac{v_{j}}{n_{j}}\textbf{X}_{i}^{j}\textbf{w}^{k}|^{s-1}\circ\texttt{sign}(\frac{v_{j}}{n_{j}}\textbf{X}_{i}^{j}\textbf{w}^{k})]$. *

*        Case 1: $0<p<1$ *

*         $\textbf{u}^{k}=|\textbf{w}^{k}|^{2-p}\circ\textbf{v}^{k},$ *

*         $\textbf{w}^{k+1}=\frac{\textbf{u}^{k}}{\parallel\textbf{u}^{k}\parallel_{p}}.$ *

*        Case 2: $p=1$ *

*         $j=\arg\max\nolimits_{i\in[1,w]}|v_{i}^{k}|,$ *

*         $w_{i}^{k+1}=\left\{\begin{aligned} \texttt{sign}(v_{j}^{k}),\ i=j,\\ 0,\ i\neq j.\\ \end{aligned}\right.$ *

*        Case 3: $1<p<\infty$ *

*         $q=p/(p-1),$ *

*         $\textbf{u}^{k}=|\textbf{v}^{k}|^{q-1}\circ\texttt{sign}(\textbf{v}^{k}),$ *

*         $\textbf{w}^{k+1}=\frac{\textbf{u}^{k}}{\parallel\textbf{u}^{k}\parallel_{p}}.$ *

*        Case 4: $p=\infty$ *

*         $\textbf{w}^{k+1}=\texttt{sign}(\textbf{v}^{k}).$ *

*        $f_{k+1}=\gamma\sum\limits^{n}_{i=1}\|\textbf{X}_{i}\textbf{w}^{k+1}\|_{s}^{s}+(1-\gamma)\sum\limits^{m}_{j=1}\sum\limits^{n_{j}}_{j=1}\|\frac{v_{j}}{n_{j}}\textbf{X}_{i}^{j}\textbf{w}^{k+1}\|_{s}^{s}.$ *

*        $\delta=|f_{k+1}-f_{k}|/|f_{k}|.$ *

*        $k\leftarrow k+1.$ *

*     end *while

*     $\textbf{W}\leftarrow[\textbf{W},\textbf{w}^{k}].$ *

*     $\textbf{D}=\texttt{diag}(\textbf{D},f_{k}).$ *

*     $\textbf{X}_{i}=\textbf{X}_{i}(\textbf{I}-\textbf{W}\textbf{W}^{T}),i=1,2,\cdots,n.$ *

*  end *for

3.3 Projection relaxation

In Section 3.2, we obtain $r$ pairs of optimal values and projection vectors: $(f_{1},\textbf{w}_{1}),\ldots,(f_{r},\textbf{w}_{r})$. Define the feature image of sample $\textbf{X}_{i}$ under W as

[TABLE]

Each column of $\textbf{P}_{i}$, $\textbf{y}_{j}=(\textbf{X}_{i}-{\bf{M}})\textbf{w}_{j}$, is called the principal component (vector).

Now we use a nearest neighbour classifier for face recognition. For a given testing sample X, compute its feature image, $\textbf{P}=(\textbf{X}-{\bf{M}})\textbf{W}$. Find out the nearest training sample $\textbf{X}_{i}$ $(1\leq i\leq n)$ whose feature image minimizes

[TABLE]

Such $\textbf{X}_{i}$ is output as the person to be recognized.

The distance, $\|(\textbf{P}_{i}-\textbf{P})\textbf{D}\|_{2}=\|(\textbf{X}_{i}-\textbf{X})\textbf{W}\textbf{D}\|_{2}$, is called relaxed distance between $\textbf{X}_{i}$ and X. Compared with originally defined distance, such as in [21], each projection axe $\textbf{w}_{j}$ is relaxed by $f_{j}$ in classification process, $j=1,\cdots,r$.

3.4 Restarted alternating direction search method

In the R2DPCA approach of face recognition, we need choose optimal $L_{s}$- and $L_{p}$-norms to maximize the recognition or classification rate. The traverse method will cost a huge amount of computational time. Instead, we present a restarted alternating direction search method of searching optimal values of $s$ and $p$; see Algorithm 3.3.

Algorithm 3.3**.**

Restarted alternating direction search method**

0:*  A finite range of $(s,p)$: $\Omega=\{(s_{i},p_{j})|i,j=1,2,\cdots,N\}$, and a positive number $\delta$. *

0:*  $s,t$. *

*  Step 1. Choose a starter $(s_{i}^{0},p_{j}^{0})\in\Omega$, randomly, and compute the recognition rate, denoted as $\varrho^{0}(i,j)$. *

*  Step 2. With $p_{j}=p_{j}^{0}$, find the maximal recognition rate in $\{(s_{i},p_{j}^{0})|$ $i=1,2,\cdots,N\}$, denoted as $\varrho^{1/2}(i,j)$, and denote the maximum point as $(s_{i}^{1},p_{j}^{0})$. *

*  Step 3. With $s_{i}=s^{1}_{i}$, find the maximal recognition rate in $\{(s_{i}^{1},p_{j})|j=1,2,\cdots,$ $N\}$,denoted as $\varrho^{1}(i,j)$, and denote the maximum point as $(s_{i}^{1},p_{j}^{1})$. *

*  Step 4. If $\varrho^{1}(i,j)=\varrho^{0}(i,j)$ and $(s_{i}^{1},p_{j}^{1})=(s_{i}^{0},p_{j}^{0})$, go to Step 5; otherwise, let $(s_{i}^{0},p_{j}^{0})=(s_{i}^{1},p_{j}^{1})$ be a new starter, and go to Step 2. *

*  Step 5. Find the maximal recognition rate $\varrho^{2}(i,j)$ in $\{(s_{i},p_{j})||s_{i}-s_{i}^{1}|\leq\delta,|p_{j}-p_{j}^{1}|\leq\delta\}\cap\Omega$, and denote the maximum point as $(s_{i}^{2},p_{j}^{2})$. *

*  Step 6. If $\varrho^{2}(i,j)\leq\varrho^{1}(i,j)$, output $s=s_{i}^{1}$ and $p=p_{j}^{1}$; otherwise, let $(s_{i}^{0},p_{j}^{0})=(s_{i}^{2},p_{j}^{2})$ be a new starter, and go to Step 2. *

If giving a enough large value $\delta$ in Algorithm 3.3, we can surely achieve the maximum value of recognition rate at optimal values $(s_{j}^{2},p_{j}^{2})\in\Omega$. The selecting process is indicated in Fig 1.

3.5 Mathematical theory of R2DPCA

R2DPCA is a generalization of G2DPCA [21]. As one of PCA-based methods, G2DPCA does not use the labels of data which possibly can impair class discrimination. To improve this, R2DPCA utilizes labels of training samples and variances within class to generate a relaxation vector, computes optimal projections maximizing the relaxes criterion, and thus enhances the class discrimination.

The working principle of R2DPCA can be clearly explained through a special case that $p=s=2$. The relaxed criterion with $L_{2}$-norm is also called generalized total scatter criterion, and has the form:

[TABLE]

where

[TABLE]

where $\textbf{W}=[\textbf{w}_{1},\cdots,\textbf{w}_{r}]$ has orthogonal columns and each column is unitary under $L_{p}$-norm, $v_{j}$ is the $j$-th element of the relaxation vector v. Here $\sum_{j=1}^{m}n_{j}=n$. Recall that ${\bf{M}}\in\mathbb{R}^{h\times w}$ is the mean of training samples. Let $\textbf{W}^{opt}=[\textbf{w}_{1}^{opt},\cdots,\textbf{w}_{r}^{opt}]$ be the optimal projection, where $\textbf{w}_{1}^{opt},\cdots,\textbf{w}_{r}^{opt}$ solve the optimal problem (21). These optimal projection axes are in fact the orthogonal eigenvectors of $\gamma\textbf{G}+(1-\gamma)\widetilde{\textbf{G}}$ corresponding to first $r$ largest eigenvalues. Since the matrix $\gamma\textbf{G}+(1-\gamma)\widetilde{\textbf{G}}$ is symmetric and positive semi-definite, $J(\textbf{W})$ is nonnegative.

R2DPCA with $s=p=2$ can also be seen as applying the relaxation idea to 2DPCA, and thus called relaxed 2DPCA. Algorithm 3.2 with $s=2,p=2$ is one method of processing the relaxed 2DPCA. Another method is applying the eigenvalue decomposition (Algorithm 3.4 ), as shown in Remark 2.2.

Algorithm 3.4**.**

relaxed 2DPCA.**

0:*  $\textbf{X}_{1},\textbf{X}_{2},\cdots,\textbf{X}_{n},\textbf{M},s=p=2,r,m,w.$ *

0:*  $\textbf{W},\textbf{D}$. *

*  $\textbf{v}=\texttt{relaxvec}(\textbf{X}_{1},\textbf{X}_{2},\cdots,\textbf{X}_{n},m,w);$ *

*  $\textbf{G}=\texttt{zeros}(w,w);$ *

*  Compute G and $\widetilde{\textbf{G}}$ defined as in (24). *

*  $[\textbf{W},\textbf{D}]=\texttt{eig}(\gamma\textbf{G}+(1-\gamma)\widetilde{\textbf{G}});$ *

*  $\textbf{W}=\textbf{W}(:,1:r);\textbf{D}=\texttt{diag}(\textbf{D}(1:r)).$ *

Remark 3.1**.**

If there is no label information or people don’t want to use it, then let $m=n$ (thus, $v_{j}=1/n,~{}n_{j}=1$) or $m=1$ (thus, $v_{j}=1,~{}n_{j}=n$), $\widetilde{\textbf{G}}$ defined in (24b) is exactly the total scatter matrix (24a) defined for the classic 2DPCA just like the case $\gamma=1$.

Now we focus on the relationships among 2DPCA, 2DPCA$L_{1}$, 2DPCA$L_{1}$-S , G2DPCA, and R2DPCA. It is obvious that 2DPCA and 2DPCA-$L_{1}$ are two special cases of G2DPCA. 2DPCA$L_{1}$-S originates from G2DPCA with $s=1$ and $p=1$ which leads to projection vector with only one nonzero element. Then the $L_{2}$-norm constraint is employed to fix this problem, resulting in 2DPCA$L_{1}$-S. On the other hand, G2DPCA with $s=1$ and $1<p<2$ behaves like 2DPCA$L_{1}$-S, since the $L_{p}$-norm constraint in G2DPCA behaves like the mixed-norm constraint in 2DPCA$L_{1}$-S. Applying the relaxation idea to 2DPCA, 2DPCA$L_{1}$ and 2DPCA$L_{1}$-S, we can get three special cases of R2DPCA. To get a better understanding of these relationships, we construct a relationship graph in Fig 2.

4 Experiments

In this section, we present numerical experiments to compare the proposed relaxed two dimensional principle component analysis (R2DPCA) by $L_{p}$-norm with state-of-the-art algorithms on face recognition. Three famous databases are utilized:

• •

Faces95 database (1440 images from 72 subjects, twenty images per subject),

• •

color Feret database (3025 images from 275 subjects, eleven images per subject),

• •

grey Feret database (1400 images from 200 subjects, seven images per subject).

All of face images are cropped and resized, and each image is of 80$\times$80 size. The numerical experiments are performed with MATLAB-R2016 on a personal computer with Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.4GHz (dual processor) and RAM 32GB.

Example 4.1**.**

In this experiment, we compare R2DPCA with 2DPCA, 2DPCA-$L_{1}$, 2DPCA$L_{1}$-S, and G2DPCA. We randomly select $10$ and $5$ images of each person from Faces95 database and color Feret face database as the training set, respectively, and the remaining as the testing set. As in [21] we set $\Omega=\{(s,p)|s=1.0:0.1:3.0,p=0.9:0.1:3.0\}$ for G2DPCA and R2DPCA. The parameter $\rho$ of 2DPCA$L_{1}$-S relates to the $\lambda$ in [11] via $\lambda=10^{-\rho}$ is tuned, consistent with [20]. The optimal $\rho$ value is selected from $[-3.0:0.1:3.0]$. Here, the relaxed parameter of criterion in (18) is set as $\gamma=0$ and the number of eigenfaces is fixed as $r=10$ (Other cases will be considered in Examples 4.2-4.3).

We repeat the whole procedure two times and output the average recognition rate. The face recognition rate (Accuracy) and corresponding optimal parameter are listed in Table 2 and Table 3, respectively. The reasonable trend of the classification accuracies according to different choices of $s$ and $p$ is presented in Fig 3 and Fig 4. These numerical results indicate that R2DPCA performs better than other four state-of-art algorithms.

Example 4.2**.**

In this experiment, we research the effect of the parameter $\gamma$ of R2DPCA on the classification accuracy. The first $10$ and $5$ images of each person are selected as the training sets of the Faces95 and color/gray Feret face databases, respectively; and $r=10$ features are selected.

The results with several representative values $\gamma$ are shown in Table 4 and Table 5. We can see that the Faces95 and color FERET databases are not sensitive to $\gamma$, and however, we can see the validity of the parameters $\gamma$ on the Gray Feret database.

Example 4.3**.**

In this experiment, we test the effect of numbers of chosen features on the classification accuracy. We randomly select $10,5$ images of each subject as training samples and the remaining as testing samples on the Faces95 and Color Feret databases, respectively. The whole procedure is repeated two times and the average accuracies are listed. Based on the optimal parameters $s,p$ of R2DPCA with $\gamma=0$ in Example 4.1, we set $s=1,~{}~{}p=2.2$ and $s=3,~{}~{}p=2.2$.

Fig 5 and Fig 6 show the classification accuracies of G2DPCA and R2DPCA with different feature numbers in the range of $[1,30]$ on the Faces95 database and Color Feret database, respectively. From these results, we can see that the classification accuracies of R2DPCA are higher and more stable than G2DPCA. When $k=1$ the classification accuracies of G2DPCA and R2DPCA are the same, which consists to the theory.

Example 4.4**.**

In this experiment, we research the influence of parameters $s$, $p$ on classification accuracies of R2DPCA with the case $\gamma=0$.

The training sets and testing sets just like Example 4.1. The classification accuracies with $10$ feature numbers, then the results are recorded. The procedure is repeated two times and then we take the average value. We use the optimal parameters of each databases on R2DPCA in Example 4.1. Here we fix $s=1$ and search the optimal parameters set from $p=[0.9:0.1:3.0]$ on the Faces95 database. We fix $s=3$ and search the optimal parameters set from $p=[0.9:0.1:3.0]$ on the Color Feret database. Similarly, we fix $p=2.2$ and search the optimal parameters set from $s=[1.0:0.1:3.0]$ on the Faces95 database. We fix $p=2.2$ and search the optimal parameters set from $s=[1.0:0.1:3.0]$ on the Color Feret database.

The results are shown in Fig 7 and Fig 8. From these Figures, we know that when $s=1,p=2.2$, the accuracy classification approach to maximum on the Faces95 database. When $s=3,p=2.2$, the accuracy classification approach to maximum on the Color Feret database. These results are consistent with the results of Experiment 4.1. And from these figures, we also can know that the accuracy classification do not have a stable variation trend with $s$ or $p$.

Example 4.5**.**

In section 3.4, we proposed a restarted alternating direction search method. Now we test this method on the Faces95 database and Color Feret database. From Example 4.4, we know the classification accuracies don’t have a stable trend with different $s$ or $p$. Traditionally, we need to traverse all the combinations of $s,p$, so that we can get the maximum solution. This way spends large time of calculations. Here, we test the restarted alternating direction search method.

$\Omega=\{(s,p)|s=1.0:0.1:3.0,p=0.9:,0.1:3.0\}$* is set as in [21]. We randomly use $10,5$ images of each subject as the training samples and the remaining images as the testing samples to do this restarted alternating direction search algorithm. We randomly start from four initial points, i.e., $s=1,s=1.6,p=0.9,p=1.2$ to find corresponding starters and then find the max classification accuracy of R2DPCA. We set a value $s$ (or $p$) to find a value $p$ (or $s$) and use $(s,p)$ as a starter. We set $\delta=0.3$ to control time. We do R2DPCA with $\gamma=0$ and obtain the max classification accuracy $0.9375,0.7158$ on two databases respectively. In order to have a more intuitive view of searching paths, we describe the results in Table 6 and Table 7.*

From section 3.4, if $\varrho^{1}(i,j)=\varrho^{0}(i,j)$ and $(s_{i}^{1},p_{j}^{1})=(s_{i}^{0},p_{j}^{0})$, we should go to Step 5 to do a restarted algorithm. This situation occurs in this experiment. For example, in table 6, we can see that if starter is $s=1,p=2.2$, we make a restarted algorithm in the next step. Because of a small positive value $\delta=0.3$, it can’t reach to the new starter correspond to the next max accuracy which can be seen from the second row in this table. If we set a enough bigger positive value $\delta$, it can search the optimal solution. Due to this property, we can use this algorithm to do a pre computation with databases which needed to be identified, so that we can have a general idea of whether these databases are applicable to the algorithms we proposed.

In each two databases, we should do Algorithm 3.2 462 times in traditional R2DPCA. But now we at most do Algorithm 3.2 153 times and even do 3 times at least. And there is a little difference between observed results and true accuracies.

Example 4.6**.**

This method is common to 2DPCA-like methods. Now we compare the results of algorithms including 2DPCA, 2DPCA-$L_{1}$, 2DPCA$L_{1}$-S with their relaxation results on Gray Feret database.

We randomly select $tr=3$ training samples from each subject and the remaining images as testing samples. Here, we also choose 10 feature numbers to save computational time. Then the nearest neighbour classifier is applied to do classification. Also, the procedure is repeated two times and take the average classification accuracies. Notice that the relaxed 2DPCA Algorithm 3.4 as shown in Remark 2.2 and the other two algorithms’ relaxed progression are similar with R2DPCA Algorithm 3.2. In order to be consistent with the previous experimental parameters, we set $\gamma=0$. The results see Table 8.

From the results, we can know that our proposed R2DPCA algorithm is also effective for other 2-D algorithms.

Example 4.7**.**

We test the accuracies of the Gray Feret database with the number of training samples in this example.

We randomly select $tr=3$ training samples from each subject and the remaining images as testing samples. Here, we also choose 10 feature numbers. Then the nearest neighbor classifier is applied to do classification. We do this process two times and take the average classification accuracies in Fig 9. Also, in order to be consistent with the previous experimental parameters, we set $\gamma=0$.

From this Figure, we known that the classification accuracies of relaxed versions are higher than those in original versions.

5 Conclusion

In this paper, we present a relaxed two dimensional principal component analysis (R2DPCA) approach for face recognition, with applying the label information of the training data. The R2DPCA is a generalization of 2DPCA, 2DPCA-$L_{1}$ and G2DPCA, and has higher generalization ability. Since utilizing the label information, the R2DPCA can be seen as a new supervised projection method, but it is totally different to the two-dimensional linear discriminant analysis (2DLDA)[35, 36].

Acknowledgments

This paper is supported in part by National Natural Science Foundation of China under grants 11771188 and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.