# Robust Structure Preserving Nonnegative Matrix Factorization for Dimensionality Reduction

Authors: Bingfeng Li; Yandong Tang; Zhi Han
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

As a linear dimensionality reduction method, nonnegative matrix factorization (NMF) has been widely used in many fields, such as machine learning and data mining. However, there are still two major drawbacks for NMF: (a) NMF can only perform semantic factorization in Euclidean space, and it fails to discover the intrinsic geometrical structure of high-dimensional data distribution. (b) NMF suffers from noisy data, which are commonly encountered in real-world applications. To address these issues, in this paper, we present a new robust structure preserving nonnegative matrix factorization (RSPNMF) framework. In RSPNMF, a local affinity graph and a distant repulsion graph are constructed to encode the geometrical information, and noisy data influence is alleviated by characterizing the data reconstruction term of NMF withl
2,1-norm instead of l
2-norm. With incorporation of the local and distant structure preservation regularization term into the robust NMF framework, our algorithm can discover a low-dimensional embedding subspace with the nature of structure preservation. RSPNMF is formulated as an optimization problem and solved by an effective iterative multiplicative update algorithm. Experimental results on some facial image datasets clustering show significant performance improvement of RSPNMF in comparison with the state-of-the-art algorithms.


---

## Body

## 1. Introduction

With the development of data acquisition technology, high-dimensional data is ubiquitous in many areas. Directly learning in original data space is prone to disturbance by noise; meanwhile, it is hard to discover the latent semantic structure [1], so low-dimensional representation learning becomes a fundamental problem in numerous data analysis algorithm [2–10]. As a common used tool for dimensionality reduction, matrix factorization has received more and more attention in recent years [2, 3, 11]. Unlike typical matrix factorization methods, such as QR and SVD, which seek for exact factorization of original matrix, nonnegative matrix factorization (NMF) approximately factorizes the original matrix into two low-dimensional nonnegative submatrixes. It is the nonnegative constraint which assures NMF can learn a sparse and part-based representation of the sample matrix. Psychological and physiological evidences have proved that NMF is consistent with how the brain interprets an object in a part-based, nonsubtractive way [12–14]. Recently, NMF has exhibited great success in a variety of applications, like text analysis [15], image processing [16, 17], and video representation [18].In order to detect the latent geometrical structure hidden in high-dimensional data space, many manifold learning algorithms have been proposed, for instance, isometric feature mapping (ISOMAP) [19], locally linear embedding (LLE) [20], Laplacian eigenmaps (LE) [21], and local tangent space alignment (LTSA) [22]. Motivated by the benefit of NMF and manifold learning, a graph regularized NMF (GNMF) method is proposed in [5]. Although GNMF displays some improvement on clustering performance, it suffer greatly from the following two problems in real-world applications: (a) GNMF only consider the ideal case where the input data is noise-free. When the input data samples is corrupted by noise, GNMF even leads to an inaccurate matrix factorization. (b) GNMF overemphasizes the affinity relationship between neighboring data samples; the repulsion relationship between data samples which is favorable for clustering is even neglected thoroughly.To address the abovementioned issues of GNMF, in this paper, we propose a new robust structure preserving nonnegative matrix factorization (RSPNMF) method. In RSPNMF, the data reconstruction term is formulated in the form of robustl
2,1-norm instead of l
2-norm, which will improve the model robustness to noisy data; meanwhile, the rotation invariance property of the model is also preserved. The success of GNMF relies on the assumption that two neighboring data samples share the same semantic label; it characterizes the compactness relationship among neighboring data samples. Inspired by this observation, we propose a new assumption; that is, data samples which are remote are generally semantic difference; and then we constrain that the low-dimensional representation of distant data samples in the original data space should also be far apart in embedding subspace. We consider that the preservation of neighboring data sample’s compactness together with distant data sample’s repulsion is helpful to describe the intrinsic geometrical structure of data distribution.The contributions of this paper are presented as follows:(1)
A novel nonnegative matrix factorization algorithm, named RSPNMF, is proposed, in which the intrinsic manifold structure is approximated with the local affinity and distant repulsion graphs. By incorporating the corresponding regularization term into robust NMF framework, our proposed algorithm can find the latent subspace of the original matrix with nature of manifold structure learning.
(2)
We formulate our approach as an optimization problem and provide an iterative multiplicative update algorithm to solve it. The convergence proof of the proposed algorithm is also provided.
(3)
Our experiments on several datasets for clustering validate our proposed algorithm. Compared with several state-of-the-art methods, the results of our RSPNMF show the better performance.The rest of the paper is organized as follows. In Section2, we briefly review the background of NMF and related works. In Section 3, our RSPNMF approach and the convergence proof of optimization scheme are presented. Some experimental results and comparison with other approach are presented in Section 4. We end this paper with some conclusions and discussion in Section 5.


## 2. Related Works

Given a nonnegative matrixX
=
x
1
,
…
,
x
N, x
i
∈
R
+
M
×
1 is a data sample. NMF attempts to find two low-dimensional nonnegative submatrixes whose product provides a best approximation to the original matrix X; that is, (1)
X
≈
U
V
T
s
.
t
.
U
∈
R
+
M
×
K
,
V
∈
R
+
N
×
K
.With the purpose of finding a compact factorization for the original data matrix, the number of basis component K usually is set less than M and N.One common data reconstruction function used to quantify the approximation of NMF is formulated as follows:(2)
O
=
X
-
U
V
T
F
2
,where ·
F denotes the Frobenius norm.Due to no-convexity of the objective function in (2) on both U and V, it is unrealistic to expect an algorithm to find its global minimum. An iterative multiplicative update algorithm which can find the local minima of (2) is proposed in [23], which are (3)
u
i
k
⟵
u
i
k
X
V
i
k
U
V
T
V
i
k
,
v
j
k
⟵
v
j
k
X
T
U
j
k
V
U
T
U
j
k
.As the update algorithm allows only additive combination among different basis components, it is believed that NMF can learn a sparse and part-based representation [4].Inspired by the work exhibited in [23], many of NMF variants have been proposed over the past years. In order to explicitly impose sparseness constraint, sparseness constrained NMF is proposed in [9, 24–28]; these works expect more sparse representation than the original NMF. To characterize different residual error model, two new distance metrics are introduced in [29, 30]. Similar to SVD presented in [31], a nonnegative matrix can also be factorized into three nonnegative components. In [32], a fast nonnegative matrix trifactorization is proposed, and then orthogonality constrained nonnegative matrix trifactorization is studied in [27]. Different from unsupervised NMF variants, a semisupervised nonnegative matrix decomposition algorithm is proposed in [9], where the label information of data samples was introduced to NMF as an additional constraint. To improve discriminating power of MMF, in [33, 34], Fisher’s criterion as an additional regularization term is incorporated into NMF framework. In spite of the success of the abovementioned NMF variants in various scenarios, all of them perform matrix factorization on Euclidean space.In order to integrate manifold learning into NMF framework, a graph regularized NMF (GNMF) algorithm is firstly proposed in [6]; then some GNMF variants are introduced in [35–38].In [39, 40], a dual regularized seminonnegative matrix trifactorization is proposed, in which data samples and feature samples manifold structure are also considered. Although the graph-based NMF algorithm incorporates the prior knowledge of manifold structure into NMF framework, overemphasizing local structure preservation sometimes will degrade its performance. A structure preserving nonnegative matrix factorization (SPNMF) is proposed in [10], in which the intrinsic manifold is approximated with local and distant graphs. Experiment results demonstrate SPNMF outperforms GNMF and its variants.All the abovementioned NMF variants are based on the premise that the input data are noise-free; it is impractical to real-world data samples which are often corrupted by noise. For boosting NMF robustness to outliers,l
2,1-norm based NMF reconstruction function is applied in [41, 42]. In [43], l
2,1-norm and l
1-norm are used to formulate the graph embedding and data reconstruction function, which endows the algorithm with robustness to unreliable graphs and noisy labels.Among the existing methods, SPNMF is closely relevant to our method and performs better on some specific datasets. There are still some differences between our proposed algorithm and SPNMF. Firstly, the Laplacian matrix of distant repulsion graph from SPNMF needs to be updated in each iteration, which is time consuming. Additionally, to learn a better parts-based representation, the basis vector of SPNMF is required to be as orthogonal as possible. Finally, SPNMF is sensitive to noisy data points.


## 3. The Proposed Method

Our goal is to perform matrix factorization on the hidden semantic subspace. To this end, we impose the local affinity and distant repulsion structure preservation on the new robust NMF framework. In this section, we will first explain the structure preservation and robust data reconstruction terms of RSPNMF, followed by its optimization algorithm. The convergence of the proposed algorithm is proved finally.


### 3.1. Problem Formulation

Despite the diverse motivations of various dimensionality reduction algorithms, most of them can be explained within a graph embedding framework. LetG
=
{
X
,
E
} be an undirected weighted graph, where vertex set X corresponds to a dataset and E
∈
R
+
N
×
N is an affinity matrix whose elements measure the similarity of each pair of vertices. In graph embedding algorithm, graph G characterizes the prior knowledge of geometric structure from data distribution.If we use thek
1-nearest neighboring graph G to characterize the local structure of data distribution, the weight matrix W
=
w
i
j
N
×
N can be defined as follows:(4)
w
i
j
=
e
-
x
i
-
x
j
2
/
σ
2
,
if
x
i
∈
N
k
1
x
j
o
r
x
j
∈
N
k
1
x
i
,
0
,
otherwise
,where σ is the bandwidth parameter and N
k
1
(
x
i
) denotes the set of k
1 nearest neighbors of x
i.Local invariance assumption which encourages neighboring data pairs in the original space to be still close in the low-dimensional embedding subspace can be formulated as(5)
m
i
n
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
,where v
i is the corresponding low-dimensional representation of x
i.By defining the diagonal matrixD and Laplacian matrix L of graph G are as (6)
L
=
D
-
W
,
d
i
i
=
∑
i
≠
j
w
i
j
,
∀
i
,we can rewrite (5) as(7)
m
i
n
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
=
m
i
n
V
∑
i
=
1
N
v
i
T
v
i
d
i
i
-
∑
i
,
j
=
1
N
v
i
T
v
j
w
i
j
=
m
i
n
V
T
r
V
T
D
V
-
T
r
V
T
W
V
=
m
i
n
V
T
r
V
T
L
V
,where T
r
(
) denotes the trace of a matrix.Equation (7) essentially holds the smoothness of dimensionality reduction process. It is based on two vital assumptions; firstly, the neighboring data samples in high-dimensional space are semantic similarity, and secondly the preservation of affinity structure plays great role for low-dimensional representation.Local invariance assumption essentially exploits the favorite relationship among similar data samples under unsupervised condition; however, it ignores unfavorite relationship between divergent data pairs. In this paper, we conjecture that the distant data pairs are always semantic difference. A new distant neighboring graphG
C
=
{
X
,
E
C
} which is used to describe the repulsion relationship between dissimilar data pairs is also constructed. The corresponding weight matrix W
C
=
w
i
j
c
N
×
N is defined as follows:(8)
w
i
j
c
=
e
-
x
i
-
x
j
2
/
σ
2
,
if
x
i
∈
D
k
2
x
j
o
r
x
j
∈
D
k
2
x
i
,
1
,
otherwise
,where D
k
2
(
x
i
) denotes the remotest k
2 data samples of x
i in givensample dataset.From (8), we can we can find that the more the distance between x
i and x
j, the smaller the value of w
i
j
c. If we want the corresponding low-dimensional representation v
i and v
j are also faraway, the value of v
i
-
v
j
2 must be as large as possible. Inspired by this observation, we propose a new distant repulsion assumption, which emphasizes that the low-dimensional counterpart of distant data pairs is also far apart in embedding subspace. This distant repulsion constraint is characterized by the objective function defined in (9)
m
a
x
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
c
.The distant repulsion assumption in (9) can maintain the repellent relationship between each remote data pair during dimensionality reduction, which can be rewritten as(10)
m
a
x
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
c
=
m
a
x
V
∑
i
=
1
N
v
i
T
v
i
d
i
i
c
-
∑
i
,
j
=
1
N
v
i
T
v
j
w
i
j
c
=
m
a
x
V
T
r
V
T
D
C
V
-
T
r
V
T
W
C
V
=
m
a
x
V
T
r
V
T
L
C
V
,where D
C and L
C are the diagonal matrix and Laplacian matrix of graph G
C, respectively, and L
C
=
D
C
-
W
C, d
i
i
C
=
∑
i
≠
j
w
i
j
C for any given i.With the assumption that semantic label of distant data pairs is certainly different, (10) actually emphasizes the unfavorite relationship among dissimilar data samples.


### 3.2. Algorithm

Although (7) and (10) together can better preserve the manifold structure of data distribution, their optimized object is opposite, one for minimum and the other for maximum. To address this issue, we divide the basis matrix U into two parts; that is,(11)
U
=
U
1
,
U
2
;correspondingly, the coefficient matrix V is also partitioned into two parts:(12)
V
=
V
1
,
V
2
.Then, the original matrix X can be reconstructed with U
1
,
V
1 and U
2
,
V
2 in an additive manner:(13)
X
≈
U
1
,
U
2
V
1
,
V
2
T
=
U
1
V
1
T
+
U
2
V
2
T
,where U
1
∈
R
+
M
×
Q, V
1
∈
R
+
N
×
Q, U
2
∈
R
+
M
×
(
K
-
Q
), V
2
∈
R
+
N
×
(
K
-
Q
), and K
≥
Q. V
1 and V
2 are the coefficient matrix of the original matrix X on U
1 and U
2, respectively.If we restrict the fact that the submatrixV
1 retains the property of graph G and avoids the property of graph G
C, then submatrix V
1 is the low-dimensional representation of matrix X with local invariance assumption. With the added constraint that the property of basic matrices U
1 and U is consistent, (7) can be translated as(14)
m
i
n
V
1
∑
i
,
j
=
1
N
v
i
1
-
v
j
1
2
W
i
j
=
m
i
n
V
1
T
r
V
1
T
L
V
1
.On the contrary, we assume that submatrix V
2 preserves the property of graph G
C and avoids the property of graph G; meanwhile, we constrain the fact that the property of matrix U
2 is complementary to matrix U. This will lead to the fact that the calculation of finding maximum value in (11) is transformed into the calculation of finding minimum value showed in (15)
m
i
n
V
2
∑
i
,
j
=
1
N
v
i
2
-
v
j
2
2
W
i
j
c
=
m
i
n
V
2
T
r
V
2
T
L
C
V
2
.In (2), the error of each data sample entered the objective function of NMF with squared residue, a few outliers with large errors will easily dominate the objective function. In this paper, we employ l
2,1-norm to define the data reconstruction term of NMF framework. Compared with l
2-norm, l
2,1-norm can keep the feature rotation invariance together with minimizing outlier impact. The l
2,1-norm of matrix X can be defined as(16)
X
2,1
=
∑
i
=
1
N
x
i
2
;correspondingly, the l
2,1 robust loss function of NMF can be formulated as(17)
O
=
X
-
U
V
T
2,1
.With the purpose of finding a robust low-dimensional representation with the nature of structure preservation, we explicitly incorporate (14) and (15) into the robust NMF framework. The objective function of the proposed RSPNMF is described as follows:(18)
O
L
=
X
-
U
V
T
2,1
2
+
α
T
r
V
1
T
L
V
1
+
β
T
r
V
2
T
L
C
V
2
s
.
t
.
U
=
U
1
,
U
2
≥
0
,
V
=
V
1
,
V
2
≥
0
,
α
≥
0
,
β
≥
0
.From (18), we can find that RSPNMF divides the ultimate matrix factorization framework into two parts, one for structure preservation and the other for data reconstruction. The constants α and β are the tradeoff parameters to balance the contribution of different regularization terms. When β
=
0, RSPNMF degenerates to GNMF [6]. When α
=
0 and β
=
0, RSPNMF even degrade to the same form of NMF [23].Because the objective function showed in (18) is not convex on both U and V together, in the following, we will manage to achieve its local minima with multiplicative iterative update algorithm.Letψ
i
k and ϕ
j
k be the Lagrange multiplier for constraints u
i
k
≥
0 and v
j
k
≥
0, respectively, and Ψ
=
ψ
i
k and Φ
=
ϕ
j
k are corresponding Lagrange multiplier matrix and the Lagrangian function L of O
L can be expressed as(19)
L
=
O
L
+
T
r
Ψ
U
T
+
T
r
Φ
V
T
.The partial derivatives ofL with respect to U and V
1, V
2 are as follows:(20)
∂
L
∂
U
=
-
2
Λ
⊙
X
V
+
2
Λ
⊙
U
V
T
V
+
ψ
,
∂
L
∂
V
1
=
-
2
Λ
⊙
X
T
U
1
+
2
Λ
⊙
U
1
V
1
T
T
U
1
+
2
Λ
⊙
U
2
V
2
T
T
U
1
+
2
α
L
V
1
+
Φ
,
∂
L
∂
V
2
=
-
2
Λ
⊙
X
T
U
2
+
2
Λ
⊙
U
1
V
1
T
T
U
2
+
2
Λ
⊙
U
2
V
2
T
T
U
2
+
2
β
L
C
V
2
+
Φ
,where Λ
=
1
/
X
-
U
V
T
·
1
I
M
×
1
,
…
,
1
/
X
-
U
V
T
·
N
I
M
×
1 and I
M
×
1 is a column vector with all elements equal to one.Using the KKT conditionsψ
i
k
u
i
k
=
0, ϕ
j
k
v
j
k
1
=
0, and ϕ
j
k
v
j
k
2
=
0, we obtain the following multiplicative update rules for RSPNMF:(21)
u
i
k
⟵
u
i
k
Λ
⊙
X
V
i
k
Λ
⊙
U
V
T
V
i
k
,
(22)
v
j
k
⟵
v
j
k
Λ
⊙
X
T
U
+
α
W
V
1
β
W
C
V
2
j
k
Λ
⊙
U
V
T
T
U
+
α
D
V
1
β
D
C
V
2
j
k
.To avoid the scale transfer problem which was pointed out in [44], we normalize U and V as follows:(23)
V
i
·
t
+
1
⟵
V
i
·
t
+
1
×
U
·
i
t
+
1
2
,
U
·
i
t
+
1
⟵
U
·
i
t
+
1
U
·
i
t
+
1
2
.The total procedure for solving the minimization problem of (18) is presented in Algorithm 1, from which we can obtain a low-dimensional representation of input data samples.Algorithm 1:RSPNMF algorithm description.
(1) InitializeU
0, V
0 as arbitrary non-negative matrices.
(2) Fort
=
0,1
,
2
,
…
,
T, Do
(1) For givenV
=
V
t, update the basis matrix U
t
+
1 according to (21).
(2) For givenU
=
U
t, update the basis matrix V
t
+
1 according to (22).
(3) NormalizeU
t
+
1 and V
t
+
1 according to (23).
(3) OutputU
=
U
T and V
=
V
T.


### 3.3. Convergence Analysis

Due to the objective function shown in (18) is nonconvex, only local optimal solutions can be obtained. In this section, we use the property of auxiliary function approach to prove the objective function O
L is nonincreasing under the update formulas in (21) and (22).Convergence Analysis for
ULemma 1.
If there exists an auxiliary functionG
u
,
u
t for F
u which satisfies G
u
,
u
t
≥
F
u and G
u
,
u
=
F
u, then F
u is nonincreasing under the update(24)
u
t
+
1
=
a
r
g
u
m
i
n
G
u
,
u
t
.Proof.
ConsiderF
u
t
+
1
≤
G
u
t
+
1
,
u
t
≤
G
u
t
,
u
t
=
F
u
t
.In order to prove the objective functionO
L is nonincreasing under the update formula in (21), we only need to construct a proper auxiliary function G
a
b
U
,
U
t which satisfies G
a
b
U
,
U
t
≥
F
a
b
U and G
a
b
U
,
U
=
F
a
b
U, where F
U denotes the part of O
L which is only relevant to U and F
a
b
U is the (
a
,
b
)th element of F
U.The Taylor series expansion ofF
a
b
U can be formulated as(25)
F
a
b
U
=
F
a
b
U
t
+
F
a
b
′
U
t
u
a
b
-
u
a
b
t
+
1
2
F
a
b
′
′
U
t
u
a
b
-
u
a
b
t
2
,where(26)
F
a
b
′
U
t
=
∂
O
L
∂
U
a
b
=
2
Λ
⊙
U
V
T
V
-
2
Λ
⊙
X
V
a
b
,
F
a
b
′
′
U
t
=
∂
F
a
b
′
U
t
∂
U
a
b
=
Λ
V
⊙
V
a
b
.Lemma 2.
Function(27)
G
a
b
U
,
U
t
=
F
a
b
U
t
+
F
a
b
′
U
t
u
a
b
-
u
a
b
t
+
Λ
⊙
U
V
T
V
a
b
u
a
b
t
u
a
b
-
u
a
b
t
2is an auxiliary function for F
a
b
U which satisfies G
a
b
U
,
U
t
≥
F
a
b
U and G
a
b
U
,
U
t
=
F
a
b
U.Proof.
Obviously,G
a
b
U
,
U
=
F
a
b
U is true; we only need to prove G
a
b
U
,
U
t
≥
F
a
b
U. For this purpose, we compare G
a
b
U
,
U
t with the Taylor series expansion of F
a
b
U.
Since(28)
Λ
V
⊙
V
a
b
u
a
b
t
=
∑
j
Λ
a
j
v
j
b
v
j
b
u
a
b
t
,
Λ
⊙
U
V
T
V
a
b
=
∑
j
Λ
a
j
v
j
b
∑
k
u
a
k
v
j
k
=
∑
j
Λ
a
j
v
j
b
u
a
b
t
v
j
b
+
∑
k
≠
b
u
a
k
v
j
k
,we can find that Λ
V
⊙
V
a
b
a
b
u
a
b
t
≤
Λ
⊙
U
V
T
V
a
b, and then we can conclude that G
a
b
U
,
U
t
≥
F
a
b
U.From Lemma2, we can find that F
a
b
U is nonincreasing under the update rule in (21).Let∂
G
a
b
U
,
U
t
/
∂
u
a
b
=
0; there must be u
∗, which makes G
a
b
U
,
U
t attain the minima, so we can obtain that (29)
u
∗
=
u
a
b
t
-
u
a
b
t
F
a
b
′
u
a
b
t
2
Λ
⊙
U
V
T
V
a
b
=
u
a
b
t
-
u
a
b
t
2
Λ
⊙
U
V
T
V
-
2
Λ
⊙
X
V
a
b
2
Λ
⊙
U
V
T
V
a
b
=
u
a
b
Λ
⊙
X
V
a
b
Λ
⊙
U
V
T
V
a
b
.Convergence Analysis for
V
1. Similar to the convergence analysis for U, the first and the second derivate of O
L with respect to V
1 can be described as(30)
F
a
b
′
V
1
=
∂
O
L
∂
V
1
a
b
=
-
2
Λ
⊙
X
T
U
1
+
2
Λ
⊙
U
V
T
T
U
1
+
2
D
-
W
V
1
a
b
,
F
a
b
′
′
V
1
=
∂
F
a
b
′
V
1
∂
V
1
a
b
=
Λ
T
U
1
⊙
U
1
a
b
+
2
D
-
W
a
a
.The Taylor series expansion ofF
a
b
V
1 is as follows:(31)
F
a
b
V
1
=
F
a
b
V
1
t
+
F
a
b
′
V
1
t
v
a
b
1
-
v
a
b
1
t
+
1
2
F
a
b
′
′
V
1
t
v
a
b
1
-
v
a
b
1
t
2
.Now, we construct a new auxiliary function forF
a
b
V
1: (32)
G
a
b
V
1
,
V
1
t
=
F
a
b
V
1
t
+
F
a
b
′
V
1
t
v
a
b
1
-
v
a
b
1
t
+
Λ
⊙
U
V
T
T
U
1
+
α
D
V
1
a
b
v
a
b
1
t
v
a
b
1
-
v
a
b
1
t
2
,and we prove which satisfies G
a
b
V
1
,
V
1
t
≥
F
a
b
V
1 and G
a
b
V
1
,
V
1
≥
F
a
b
V
1
.Proof.
Obviously,G
a
b
V
1
,
V
1
=
F
a
b
V
1 is true; we only need to prove G
a
b
V
1
,
V
1
t
≥
F
a
b
V
1.
Since(33)
Λ
T
U
1
⊙
U
1
a
b
v
a
b
1
t
≤
Λ
⊙
U
V
T
T
U
1
a
b
,
D
-
W
a
a
v
a
b
1
t
≤
D
V
1
a
b
,we can deduce that G
v
a
b
1
,
v
a
b
1
t
≥
F
a
b
v
a
b
1
.
Let∂
G
a
b
V
1
,
V
1
t
/
∂
v
a
b
1
=
0; we can find that(34)
F
a
b
V
1
t
+
2
Λ
⊙
U
V
T
T
U
1
+
α
D
V
a
b
1
a
b
v
a
b
1
t
v
a
b
1
-
v
a
b
1
t
=
0
.Finally, we have (35)
v
a
b
1
⟵
v
a
b
1
Λ
⊙
X
T
U
1
+
α
W
V
1
a
b
Λ
⊙
U
V
T
T
U
1
+
α
D
V
1
a
b
.The convergence analysis forV
2 is similar to the convergence analysis for V
1, and it is omitted.Computational Complexity Analysis. With the update rules in (21) and (22), it is not hard to calculate the extra computational cost of our proposed algorithm in each iteration, and we summarize the result in Table 1. It is important to note that matrix W is sparse, because the local graph G is constructed with k
1 neighbors. The average nonzero element on each row of W is k
1. Thus, we only need N
k
1
Q and flam (a floating point addition and multiplication) to compute W
V
1.Table 1
Computational operation counts for updatingUand V in each iteration.
Addition
Multiplication
Division
Overall
UpdatingU (21)
3
M
N
K
-
M
(
N
+
2
K
)
3
M
N
K
+
2
M
N
M
K
O
(
M
N
K
)
UpdatingV (22)
2
N
K
(
M
-
1
)
+
M
N
(
K
-
1
)
+
N
Q
(
k
1
-
1
)
+
N
2
(
K
-
Q
)
3
M
N
K
+
N
Q
(
k
1
+
1
)
+
N
2
(
K
-
Q
)
N
K
O
(
M
N
K
)
NormalizingU and V
(
M
-
1
)
N
M
N
M
K
O
(
M
N
)
Overall
O
(
M
N
K
)
O
(
M
N
K
)
O
(
(
M
+
N
)
K
)
Addition: a floating-point addition; multiplication: a floating-point multiplication; division: a floating-point division.Besides the multiplicative updates, RSPNMF also needsO
(
M
N
2
) to construct the local affinity graph G and distant repulsion graph G
C. If the multiplicative updates stop after t iterations, the overall cost for RSPNMF is O
(
t
M
N
K
+
2
M
N
2
). The overall cost for GNMF is O
(
t
M
N
K
+
M
N
2
) and the overall cost for SPNMF is O
(
t
(
M
+
N
)
N
K
+
M
N
2
).


## 3.1. Problem Formulation

Despite the diverse motivations of various dimensionality reduction algorithms, most of them can be explained within a graph embedding framework. LetG
=
{
X
,
E
} be an undirected weighted graph, where vertex set X corresponds to a dataset and E
∈
R
+
N
×
N is an affinity matrix whose elements measure the similarity of each pair of vertices. In graph embedding algorithm, graph G characterizes the prior knowledge of geometric structure from data distribution.If we use thek
1-nearest neighboring graph G to characterize the local structure of data distribution, the weight matrix W
=
w
i
j
N
×
N can be defined as follows:(4)
w
i
j
=
e
-
x
i
-
x
j
2
/
σ
2
,
if
x
i
∈
N
k
1
x
j
o
r
x
j
∈
N
k
1
x
i
,
0
,
otherwise
,where σ is the bandwidth parameter and N
k
1
(
x
i
) denotes the set of k
1 nearest neighbors of x
i.Local invariance assumption which encourages neighboring data pairs in the original space to be still close in the low-dimensional embedding subspace can be formulated as(5)
m
i
n
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
,where v
i is the corresponding low-dimensional representation of x
i.By defining the diagonal matrixD and Laplacian matrix L of graph G are as (6)
L
=
D
-
W
,
d
i
i
=
∑
i
≠
j
w
i
j
,
∀
i
,we can rewrite (5) as(7)
m
i
n
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
=
m
i
n
V
∑
i
=
1
N
v
i
T
v
i
d
i
i
-
∑
i
,
j
=
1
N
v
i
T
v
j
w
i
j
=
m
i
n
V
T
r
V
T
D
V
-
T
r
V
T
W
V
=
m
i
n
V
T
r
V
T
L
V
,where T
r
(
) denotes the trace of a matrix.Equation (7) essentially holds the smoothness of dimensionality reduction process. It is based on two vital assumptions; firstly, the neighboring data samples in high-dimensional space are semantic similarity, and secondly the preservation of affinity structure plays great role for low-dimensional representation.Local invariance assumption essentially exploits the favorite relationship among similar data samples under unsupervised condition; however, it ignores unfavorite relationship between divergent data pairs. In this paper, we conjecture that the distant data pairs are always semantic difference. A new distant neighboring graphG
C
=
{
X
,
E
C
} which is used to describe the repulsion relationship between dissimilar data pairs is also constructed. The corresponding weight matrix W
C
=
w
i
j
c
N
×
N is defined as follows:(8)
w
i
j
c
=
e
-
x
i
-
x
j
2
/
σ
2
,
if
x
i
∈
D
k
2
x
j
o
r
x
j
∈
D
k
2
x
i
,
1
,
otherwise
,where D
k
2
(
x
i
) denotes the remotest k
2 data samples of x
i in givensample dataset.From (8), we can we can find that the more the distance between x
i and x
j, the smaller the value of w
i
j
c. If we want the corresponding low-dimensional representation v
i and v
j are also faraway, the value of v
i
-
v
j
2 must be as large as possible. Inspired by this observation, we propose a new distant repulsion assumption, which emphasizes that the low-dimensional counterpart of distant data pairs is also far apart in embedding subspace. This distant repulsion constraint is characterized by the objective function defined in (9)
m
a
x
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
c
.The distant repulsion assumption in (9) can maintain the repellent relationship between each remote data pair during dimensionality reduction, which can be rewritten as(10)
m
a
x
V
1
2
∑
i
,
j
=
1
N
v
i
-
v
j
2
w
i
j
c
=
m
a
x
V
∑
i
=
1
N
v
i
T
v
i
d
i
i
c
-
∑
i
,
j
=
1
N
v
i
T
v
j
w
i
j
c
=
m
a
x
V
T
r
V
T
D
C
V
-
T
r
V
T
W
C
V
=
m
a
x
V
T
r
V
T
L
C
V
,where D
C and L
C are the diagonal matrix and Laplacian matrix of graph G
C, respectively, and L
C
=
D
C
-
W
C, d
i
i
C
=
∑
i
≠
j
w
i
j
C for any given i.With the assumption that semantic label of distant data pairs is certainly different, (10) actually emphasizes the unfavorite relationship among dissimilar data samples.


## 3.2. Algorithm

Although (7) and (10) together can better preserve the manifold structure of data distribution, their optimized object is opposite, one for minimum and the other for maximum. To address this issue, we divide the basis matrix U into two parts; that is,(11)
U
=
U
1
,
U
2
;correspondingly, the coefficient matrix V is also partitioned into two parts:(12)
V
=
V
1
,
V
2
.Then, the original matrix X can be reconstructed with U
1
,
V
1 and U
2
,
V
2 in an additive manner:(13)
X
≈
U
1
,
U
2
V
1
,
V
2
T
=
U
1
V
1
T
+
U
2
V
2
T
,where U
1
∈
R
+
M
×
Q, V
1
∈
R
+
N
×
Q, U
2
∈
R
+
M
×
(
K
-
Q
), V
2
∈
R
+
N
×
(
K
-
Q
), and K
≥
Q. V
1 and V
2 are the coefficient matrix of the original matrix X on U
1 and U
2, respectively.If we restrict the fact that the submatrixV
1 retains the property of graph G and avoids the property of graph G
C, then submatrix V
1 is the low-dimensional representation of matrix X with local invariance assumption. With the added constraint that the property of basic matrices U
1 and U is consistent, (7) can be translated as(14)
m
i
n
V
1
∑
i
,
j
=
1
N
v
i
1
-
v
j
1
2
W
i
j
=
m
i
n
V
1
T
r
V
1
T
L
V
1
.On the contrary, we assume that submatrix V
2 preserves the property of graph G
C and avoids the property of graph G; meanwhile, we constrain the fact that the property of matrix U
2 is complementary to matrix U. This will lead to the fact that the calculation of finding maximum value in (11) is transformed into the calculation of finding minimum value showed in (15)
m
i
n
V
2
∑
i
,
j
=
1
N
v
i
2
-
v
j
2
2
W
i
j
c
=
m
i
n
V
2
T
r
V
2
T
L
C
V
2
.In (2), the error of each data sample entered the objective function of NMF with squared residue, a few outliers with large errors will easily dominate the objective function. In this paper, we employ l
2,1-norm to define the data reconstruction term of NMF framework. Compared with l
2-norm, l
2,1-norm can keep the feature rotation invariance together with minimizing outlier impact. The l
2,1-norm of matrix X can be defined as(16)
X
2,1
=
∑
i
=
1
N
x
i
2
;correspondingly, the l
2,1 robust loss function of NMF can be formulated as(17)
O
=
X
-
U
V
T
2,1
.With the purpose of finding a robust low-dimensional representation with the nature of structure preservation, we explicitly incorporate (14) and (15) into the robust NMF framework. The objective function of the proposed RSPNMF is described as follows:(18)
O
L
=
X
-
U
V
T
2,1
2
+
α
T
r
V
1
T
L
V
1
+
β
T
r
V
2
T
L
C
V
2
s
.
t
.
U
=
U
1
,
U
2
≥
0
,
V
=
V
1
,
V
2
≥
0
,
α
≥
0
,
β
≥
0
.From (18), we can find that RSPNMF divides the ultimate matrix factorization framework into two parts, one for structure preservation and the other for data reconstruction. The constants α and β are the tradeoff parameters to balance the contribution of different regularization terms. When β
=
0, RSPNMF degenerates to GNMF [6]. When α
=
0 and β
=
0, RSPNMF even degrade to the same form of NMF [23].Because the objective function showed in (18) is not convex on both U and V together, in the following, we will manage to achieve its local minima with multiplicative iterative update algorithm.Letψ
i
k and ϕ
j
k be the Lagrange multiplier for constraints u
i
k
≥
0 and v
j
k
≥
0, respectively, and Ψ
=
ψ
i
k and Φ
=
ϕ
j
k are corresponding Lagrange multiplier matrix and the Lagrangian function L of O
L can be expressed as(19)
L
=
O
L
+
T
r
Ψ
U
T
+
T
r
Φ
V
T
.The partial derivatives ofL with respect to U and V
1, V
2 are as follows:(20)
∂
L
∂
U
=
-
2
Λ
⊙
X
V
+
2
Λ
⊙
U
V
T
V
+
ψ
,
∂
L
∂
V
1
=
-
2
Λ
⊙
X
T
U
1
+
2
Λ
⊙
U
1
V
1
T
T
U
1
+
2
Λ
⊙
U
2
V
2
T
T
U
1
+
2
α
L
V
1
+
Φ
,
∂
L
∂
V
2
=
-
2
Λ
⊙
X
T
U
2
+
2
Λ
⊙
U
1
V
1
T
T
U
2
+
2
Λ
⊙
U
2
V
2
T
T
U
2
+
2
β
L
C
V
2
+
Φ
,where Λ
=
1
/
X
-
U
V
T
·
1
I
M
×
1
,
…
,
1
/
X
-
U
V
T
·
N
I
M
×
1 and I
M
×
1 is a column vector with all elements equal to one.Using the KKT conditionsψ
i
k
u
i
k
=
0, ϕ
j
k
v
j
k
1
=
0, and ϕ
j
k
v
j
k
2
=
0, we obtain the following multiplicative update rules for RSPNMF:(21)
u
i
k
⟵
u
i
k
Λ
⊙
X
V
i
k
Λ
⊙
U
V
T
V
i
k
,
(22)
v
j
k
⟵
v
j
k
Λ
⊙
X
T
U
+
α
W
V
1
β
W
C
V
2
j
k
Λ
⊙
U
V
T
T
U
+
α
D
V
1
β
D
C
V
2
j
k
.To avoid the scale transfer problem which was pointed out in [44], we normalize U and V as follows:(23)
V
i
·
t
+
1
⟵
V
i
·
t
+
1
×
U
·
i
t
+
1
2
,
U
·
i
t
+
1
⟵
U
·
i
t
+
1
U
·
i
t
+
1
2
.The total procedure for solving the minimization problem of (18) is presented in Algorithm 1, from which we can obtain a low-dimensional representation of input data samples.Algorithm 1:RSPNMF algorithm description.
(1) InitializeU
0, V
0 as arbitrary non-negative matrices.
(2) Fort
=
0,1
,
2
,
…
,
T, Do
(1) For givenV
=
V
t, update the basis matrix U
t
+
1 according to (21).
(2) For givenU
=
U
t, update the basis matrix V
t
+
1 according to (22).
(3) NormalizeU
t
+
1 and V
t
+
1 according to (23).
(3) OutputU
=
U
T and V
=
V
T.


## 3.3. Convergence Analysis

Due to the objective function shown in (18) is nonconvex, only local optimal solutions can be obtained. In this section, we use the property of auxiliary function approach to prove the objective function O
L is nonincreasing under the update formulas in (21) and (22).Convergence Analysis for
ULemma 1.
If there exists an auxiliary functionG
u
,
u
t for F
u which satisfies G
u
,
u
t
≥
F
u and G
u
,
u
=
F
u, then F
u is nonincreasing under the update(24)
u
t
+
1
=
a
r
g
u
m
i
n
G
u
,
u
t
.Proof.
ConsiderF
u
t
+
1
≤
G
u
t
+
1
,
u
t
≤
G
u
t
,
u
t
=
F
u
t
.In order to prove the objective functionO
L is nonincreasing under the update formula in (21), we only need to construct a proper auxiliary function G
a
b
U
,
U
t which satisfies G
a
b
U
,
U
t
≥
F
a
b
U and G
a
b
U
,
U
=
F
a
b
U, where F
U denotes the part of O
L which is only relevant to U and F
a
b
U is the (
a
,
b
)th element of F
U.The Taylor series expansion ofF
a
b
U can be formulated as(25)
F
a
b
U
=
F
a
b
U
t
+
F
a
b
′
U
t
u
a
b
-
u
a
b
t
+
1
2
F
a
b
′
′
U
t
u
a
b
-
u
a
b
t
2
,where(26)
F
a
b
′
U
t
=
∂
O
L
∂
U
a
b
=
2
Λ
⊙
U
V
T
V
-
2
Λ
⊙
X
V
a
b
,
F
a
b
′
′
U
t
=
∂
F
a
b
′
U
t
∂
U
a
b
=
Λ
V
⊙
V
a
b
.Lemma 2.
Function(27)
G
a
b
U
,
U
t
=
F
a
b
U
t
+
F
a
b
′
U
t
u
a
b
-
u
a
b
t
+
Λ
⊙
U
V
T
V
a
b
u
a
b
t
u
a
b
-
u
a
b
t
2is an auxiliary function for F
a
b
U which satisfies G
a
b
U
,
U
t
≥
F
a
b
U and G
a
b
U
,
U
t
=
F
a
b
U.Proof.
Obviously,G
a
b
U
,
U
=
F
a
b
U is true; we only need to prove G
a
b
U
,
U
t
≥
F
a
b
U. For this purpose, we compare G
a
b
U
,
U
t with the Taylor series expansion of F
a
b
U.
Since(28)
Λ
V
⊙
V
a
b
u
a
b
t
=
∑
j
Λ
a
j
v
j
b
v
j
b
u
a
b
t
,
Λ
⊙
U
V
T
V
a
b
=
∑
j
Λ
a
j
v
j
b
∑
k
u
a
k
v
j
k
=
∑
j
Λ
a
j
v
j
b
u
a
b
t
v
j
b
+
∑
k
≠
b
u
a
k
v
j
k
,we can find that Λ
V
⊙
V
a
b
a
b
u
a
b
t
≤
Λ
⊙
U
V
T
V
a
b, and then we can conclude that G
a
b
U
,
U
t
≥
F
a
b
U.From Lemma2, we can find that F
a
b
U is nonincreasing under the update rule in (21).Let∂
G
a
b
U
,
U
t
/
∂
u
a
b
=
0; there must be u
∗, which makes G
a
b
U
,
U
t attain the minima, so we can obtain that (29)
u
∗
=
u
a
b
t
-
u
a
b
t
F
a
b
′
u
a
b
t
2
Λ
⊙
U
V
T
V
a
b
=
u
a
b
t
-
u
a
b
t
2
Λ
⊙
U
V
T
V
-
2
Λ
⊙
X
V
a
b
2
Λ
⊙
U
V
T
V
a
b
=
u
a
b
Λ
⊙
X
V
a
b
Λ
⊙
U
V
T
V
a
b
.Convergence Analysis for
V
1. Similar to the convergence analysis for U, the first and the second derivate of O
L with respect to V
1 can be described as(30)
F
a
b
′
V
1
=
∂
O
L
∂
V
1
a
b
=
-
2
Λ
⊙
X
T
U
1
+
2
Λ
⊙
U
V
T
T
U
1
+
2
D
-
W
V
1
a
b
,
F
a
b
′
′
V
1
=
∂
F
a
b
′
V
1
∂
V
1
a
b
=
Λ
T
U
1
⊙
U
1
a
b
+
2
D
-
W
a
a
.The Taylor series expansion ofF
a
b
V
1 is as follows:(31)
F
a
b
V
1
=
F
a
b
V
1
t
+
F
a
b
′
V
1
t
v
a
b
1
-
v
a
b
1
t
+
1
2
F
a
b
′
′
V
1
t
v
a
b
1
-
v
a
b
1
t
2
.Now, we construct a new auxiliary function forF
a
b
V
1: (32)
G
a
b
V
1
,
V
1
t
=
F
a
b
V
1
t
+
F
a
b
′
V
1
t
v
a
b
1
-
v
a
b
1
t
+
Λ
⊙
U
V
T
T
U
1
+
α
D
V
1
a
b
v
a
b
1
t
v
a
b
1
-
v
a
b
1
t
2
,and we prove which satisfies G
a
b
V
1
,
V
1
t
≥
F
a
b
V
1 and G
a
b
V
1
,
V
1
≥
F
a
b
V
1
.Proof.
Obviously,G
a
b
V
1
,
V
1
=
F
a
b
V
1 is true; we only need to prove G
a
b
V
1
,
V
1
t
≥
F
a
b
V
1.
Since(33)
Λ
T
U
1
⊙
U
1
a
b
v
a
b
1
t
≤
Λ
⊙
U
V
T
T
U
1
a
b
,
D
-
W
a
a
v
a
b
1
t
≤
D
V
1
a
b
,we can deduce that G
v
a
b
1
,
v
a
b
1
t
≥
F
a
b
v
a
b
1
.
Let∂
G
a
b
V
1
,
V
1
t
/
∂
v
a
b
1
=
0; we can find that(34)
F
a
b
V
1
t
+
2
Λ
⊙
U
V
T
T
U
1
+
α
D
V
a
b
1
a
b
v
a
b
1
t
v
a
b
1
-
v
a
b
1
t
=
0
.Finally, we have (35)
v
a
b
1
⟵
v
a
b
1
Λ
⊙
X
T
U
1
+
α
W
V
1
a
b
Λ
⊙
U
V
T
T
U
1
+
α
D
V
1
a
b
.The convergence analysis forV
2 is similar to the convergence analysis for V
1, and it is omitted.Computational Complexity Analysis. With the update rules in (21) and (22), it is not hard to calculate the extra computational cost of our proposed algorithm in each iteration, and we summarize the result in Table 1. It is important to note that matrix W is sparse, because the local graph G is constructed with k
1 neighbors. The average nonzero element on each row of W is k
1. Thus, we only need N
k
1
Q and flam (a floating point addition and multiplication) to compute W
V
1.Table 1
Computational operation counts for updatingUand V in each iteration.
Addition
Multiplication
Division
Overall
UpdatingU (21)
3
M
N
K
-
M
(
N
+
2
K
)
3
M
N
K
+
2
M
N
M
K
O
(
M
N
K
)
UpdatingV (22)
2
N
K
(
M
-
1
)
+
M
N
(
K
-
1
)
+
N
Q
(
k
1
-
1
)
+
N
2
(
K
-
Q
)
3
M
N
K
+
N
Q
(
k
1
+
1
)
+
N
2
(
K
-
Q
)
N
K
O
(
M
N
K
)
NormalizingU and V
(
M
-
1
)
N
M
N
M
K
O
(
M
N
)
Overall
O
(
M
N
K
)
O
(
M
N
K
)
O
(
(
M
+
N
)
K
)
Addition: a floating-point addition; multiplication: a floating-point multiplication; division: a floating-point division.Besides the multiplicative updates, RSPNMF also needsO
(
M
N
2
) to construct the local affinity graph G and distant repulsion graph G
C. If the multiplicative updates stop after t iterations, the overall cost for RSPNMF is O
(
t
M
N
K
+
2
M
N
2
). The overall cost for GNMF is O
(
t
M
N
K
+
M
N
2
) and the overall cost for SPNMF is O
(
t
(
M
+
N
)
N
K
+
M
N
2
).


## 4. Experimental Results

In this section, we carry out several experiments to examine the clustering performance of our approach in comparison with some state-of-the-art methods.


### 4.1. Datasets Description

Three image datasets which are widely employed in dimensionality reduction literatures are used in our experiments.(
1) ORL Dataset. The ORL dataset contains 400 face images of 40 subjects under different lighting, facial expressions (open/closed eyes, smiling/not smiling), and facial details (with glasses or no glasses). All these images were taken against a dark homogeneous background with the subjects in an upright, frontal position. In our experiments, original images are firstly normalized in scale and orientation, such that two eyes are aligned at the same position. Then, the facial areas are cropped into the final images for clustering. Finally, each image is normalized to array of 32 × 32 pixels and reshaped to a long vector with 256 gray levels per pixel.(
2) YALE Dataset. The YALE dataset includes 165 face images collected from 15 volunteers. Each image has different facial expression and facial details. The same preprocessing as for the ORL dataset is implemented on YALE dataset.(
3) FERET Dataset. For the FERET dataset, we use 70 subjects with 6 facial images for each subject. Each image is also represented by a 1024-dimensional vector in image space.The statistic property of these datasets is listed in Table2, which includes the number of samples, the dimension of features, and the number of object classes.Table 2
Statistic of real-world datasets.
Datasets
# samples
# features
# class
ORL
400
1024
40
YALE
165
1024
15
FERET
420
1024
70As showed in Figure1, to evaluate the robustness of mentioned algorithms to noisy data, all the images were occluded with local and global noise. The local noise consists of occlusion patch of 12 × 12 pixels with random black and white dots. The global noise is replaced by Gaussian noise with mean 0 and variance 36.Figure 1
Sample images from ORL (a), YALE (b), and FERET (c).
(a)
(b)
(c)


### 4.2. Compared Methods

To investigate the clustering performance of our proposed approach, we compared it with the following methods:(i)
CanonicalK-means (Kmeans in short) clustering method as a baseline.
(ii)
K-means clustering in the principle component subspace (PCA in short).
(iii)
Nonnegative matrix factorization based clustering (NMF in short). We use theF-norm formulation.
(iv)
Graph regularized nonnegative matrix factorization based clustering withF-norm formulation (GNMF in short).
(v)
Structure preserving nonnegative matrix factorization based clustering (SPNMF in short).Except theKmeans, for each image, a new low-dimensional representation can be derived from the matrix factorization, and then we perform clustering by using K-means in the new embedding subspace.


### 4.3. Evaluation Metrics

To evaluate the clustering performance of all involved methods, we adopted two popular metrics, the accuracy (AC) and the normalized mutual information (NMI), as used in [6].The accuracy indicates one-to-one relationship between obtained cluster label and ground truth label. The accuracy is defined as follows:(36)
A
C
=
∑
i
=
1
n
δ
m
a
p
r
i
,
l
i
N
,where δ
m
a
p
r
i
,
l
i is the delta function which judges whether the ground truth l
i is equal to the obtained label r
i and m
a
p
r
i is a permutation mapping function that maps each obtained cluster label r
i to the equivalent label from the dataset.The second measure metrics is the normalized mutual information (NMI), which is used for determining the quality of clustering.LetC and C
′ denote the sets of cluster from the ground truth and the clustering method, respectively. The mutual information metric M
I
(
C
,
C
′
) is defined as follows:(37)
M
I
C
,
C
′
=
∑
c
i
∈
C
,
c
j
′
∈
C
′
p
c
j
,
c
j
′
·
l
o
g
2
p
c
j
,
c
j
′
p
c
j
·
p
c
j
′
,where p
c
j and p
c
j
′ denote the extent of a sample belonging to the clusters c
i and c
j
′, respectively. p
c
j
,
c
j
′ is a joint probability that an arbitrarily selected sample belongs to c
i as well as c
j
′ at the same time. The NMI can be defined as follows:(38)
N
M
I
C
,
C
′
=
M
I
C
,
C
′
max
⁡
H
C
,
H
C
′
,where H
C and H
C
′ are the entropies of C and C
′, respectively. It can be seen that N
M
I
C
,
C
′ ranges from 0 to 1. The larger the N
M
I is, the better the clustering result will be.


### 4.4. Clustering Results and Analysis

In order to fair examine the clustering performance of all mentioned methods, 20 test runs are conducted on the optimal parameter setting, and then the means and standard deviations are recorded. In each test,Kmeans is applied 20 times with different initialization and the best result is reported. All clustering performances on the three noise contaminated datasets are reported in Tables 3 and 4.Table 3
Clustering performance on local noise contaminated datasets. The best performances are highlighted in bold font.
(a)
Dataset
Accuracy (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
23.10 ± 1.14
26.12 ± 1.27
24.20 ± 1.03
25.05 ± 0.74
27.45 ± 5.22
33.06 ± 1.25
ORL
35.71 ± 2.74
46.99 ± 2.69
44.56 ± 1.59
37.78 ± 1.37
44.78 ± 1.7
54.01 ± 3.3
Yale
32.12 ± 2.38
37.43 ± 2.27
38.18 ± 2.18
32.45 ± 1.16
36.36 ± 1.24
44.42 ± 1.08
(b)
Dataset
Normalized mutual information (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
49.2 ± 2.56
56.85 ± 1.54
54.33 ± 1.12
50.06 ± 2.51
58.26 ± 0.89
67.44 ± 1.73
ORL
53.78 ± 2.56
65.79 ± 1.54
64.57 ± 1.12
58.65 ± 2.51
64.26 ± 0.89
71.50 ± 1.73
Yale
35.73 ± 2.49
42.72 ± 2.30
43.77 ± 1.58
38.84 ± 0.98
42.72 ± 1.12
47.56 ± 1.19Table 4
Clustering performance on global noise contaminated datasets. The best performances are highlighted in bold font.
(a)
Dataset
Accuracy (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
23.76 ± 1.68
27.14 ± 1.57
26.9 ± 1.23
25.76 ± 0.58
28.24 ± 0.73
35.30 ± 0.81
ORL
42.36 ± 3.24
52.78 ± 2.48
52.76 ± 2.85
52.85 ± 0.79
53.91 ± 1.78
61.5 ± 1.92
Yale
36.79 ± 2.24
41.21 ± 2.98
42.03 ± 2.83
40.06 ± 1.43
43.64 ± 1.39
50.30 ± 1.57
(b)
Dataset
Normalized mutual information (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
50.06 ± 2.74
58.95 ± 1.02
57.94 ± 1.28
60.33 ± 0.60
62.01 ± 0.57
68.65 ± 0.43
ORL
60.14 ± 2.77
70.33 ± 1.39
67.37 ± 1.94
66.7 ± 0.80
70.99 ± 0.74
76.97 ± 0.92
Yale
41.68 ± 2.27
44.78 ± 2.09
46.94 ± 2.37
45.94 ± 0.88
47.68 ± 1.08
54.89 ± 2.08From these results, we can derive the following observations: (a) SPNMF outperformsKmeans, PCA, and NMF in most cases. This demonstrates the importance of structure preservation in discovering the hidden factors. In addition, the performance of SPNMF is superior to GNMF. This reveals that, by leveraging the ratio of local affinity and distant repulsion regularization, SPNMF can learn a better low-dimensional representation. (b) The performance of RSPNMF exceeds SPNMF; this demonstrates the robustness of l
2,1-norm in handling noise data samples.


### 4.5. Parameter Sensitiveness

Our RSPNMF has five essential parameters: the number of nearest neighborsk
1, the number of distant neighbors k
2, the regularization parameter α, the regularization parameter β, and the number of basis vectors r. For all mentioned dimensionality reduction algorithms, the number of basis vectors r is equal to the number of clusters. The complement space and the original space reduced dimension d is set to r
/
2. For all graph-based methods, the number of nearest neighbors is set to {
3,5
,
…
,
19,21
}. For SPNMF and RSPNMF, the number of distant points is also set to {
3,5
,
…
,
19,21
}. For GNMF, the regularization parameter is set by the grid {
10
-
3
,
10
-
2
,
…
,
10
3
}. For SPNMF and RSPNMF, we tune α and β by searching the grid {
10
-
3
,
10
-
2
,
…
,
10
3
}. For SPNMF, γ is set to 1.The result of parameter sensitiveness on different datasets with local and global noise is showed in Figure2. In order to compare with GNMF, we fix the parameter β and tune the parameter α. From Figure 2, we can find that the performance of RSPNMF and SPNMF is superior to GNMF; meanwhile, they are all stable with respect to the change of parameter α. This demonstrates the importance of the distant repulsion preservation on characterizing the intrinsic distribution of data samples. Besides, RSPNMF still outperforms SPNMF when α range within {
0.001,10
}. This is because the robust data reconstruction term of RSPNMF alleviates the noisy data impact on clustering performance.Figure 2
The performance of involved methods with respect to the parameterα on different dataset. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALEThe relevant samples selection plays an important role in all graph regularized NMF method. For GNMF, which only usek
1-nearest neighboring graph to capture the local geometric structure of data distribution; meanwhile, assume that the neighboring samples are with the same semantic label. For noise corrupted data samples, this assumption may not hold. In particular, as k
1 increases, clustering performance of GNMF is more likely to fail as shown in Figure 3. For RSPNMF and SPNMF, the local and distant data samples relationship is also considered. This will definitely reinforce the algorithm robustness to k
1-nearest neighboring samples increase. In addition, from Figure 3 we can find that the introduction of l
2,1-norm improves RSPNMF clustering results.Figure 3
The performance of involved methods with respect to parameterk
1. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALE


### 4.6. Convergence and Computational Cost

The effectiveness of the proposed RSPNMF has been demonstrated by the above clustering experiments. Here, we evaluate the efficiency of RSPNMF, that is, how fast the update rules in (21) and (22) are converged.Figure4 shows the convergence curves of RSPNMF on all the three datasets. For each figure, the y-axis is the value of objective function and the x-axis denotes the iteration number. From the convergence curves, it can be observed that the convergence rate of RSPNMF is fast, and the runtime cost of RSPNMF is acceptable.Figure 4
The convergence curve of RSPNMF with 500 iterations. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALE


## 4.1. Datasets Description

Three image datasets which are widely employed in dimensionality reduction literatures are used in our experiments.(
1) ORL Dataset. The ORL dataset contains 400 face images of 40 subjects under different lighting, facial expressions (open/closed eyes, smiling/not smiling), and facial details (with glasses or no glasses). All these images were taken against a dark homogeneous background with the subjects in an upright, frontal position. In our experiments, original images are firstly normalized in scale and orientation, such that two eyes are aligned at the same position. Then, the facial areas are cropped into the final images for clustering. Finally, each image is normalized to array of 32 × 32 pixels and reshaped to a long vector with 256 gray levels per pixel.(
2) YALE Dataset. The YALE dataset includes 165 face images collected from 15 volunteers. Each image has different facial expression and facial details. The same preprocessing as for the ORL dataset is implemented on YALE dataset.(
3) FERET Dataset. For the FERET dataset, we use 70 subjects with 6 facial images for each subject. Each image is also represented by a 1024-dimensional vector in image space.The statistic property of these datasets is listed in Table2, which includes the number of samples, the dimension of features, and the number of object classes.Table 2
Statistic of real-world datasets.
Datasets
# samples
# features
# class
ORL
400
1024
40
YALE
165
1024
15
FERET
420
1024
70As showed in Figure1, to evaluate the robustness of mentioned algorithms to noisy data, all the images were occluded with local and global noise. The local noise consists of occlusion patch of 12 × 12 pixels with random black and white dots. The global noise is replaced by Gaussian noise with mean 0 and variance 36.Figure 1
Sample images from ORL (a), YALE (b), and FERET (c).
(a)
(b)
(c)


## 4.2. Compared Methods

To investigate the clustering performance of our proposed approach, we compared it with the following methods:(i)
CanonicalK-means (Kmeans in short) clustering method as a baseline.
(ii)
K-means clustering in the principle component subspace (PCA in short).
(iii)
Nonnegative matrix factorization based clustering (NMF in short). We use theF-norm formulation.
(iv)
Graph regularized nonnegative matrix factorization based clustering withF-norm formulation (GNMF in short).
(v)
Structure preserving nonnegative matrix factorization based clustering (SPNMF in short).Except theKmeans, for each image, a new low-dimensional representation can be derived from the matrix factorization, and then we perform clustering by using K-means in the new embedding subspace.


## 4.3. Evaluation Metrics

To evaluate the clustering performance of all involved methods, we adopted two popular metrics, the accuracy (AC) and the normalized mutual information (NMI), as used in [6].The accuracy indicates one-to-one relationship between obtained cluster label and ground truth label. The accuracy is defined as follows:(36)
A
C
=
∑
i
=
1
n
δ
m
a
p
r
i
,
l
i
N
,where δ
m
a
p
r
i
,
l
i is the delta function which judges whether the ground truth l
i is equal to the obtained label r
i and m
a
p
r
i is a permutation mapping function that maps each obtained cluster label r
i to the equivalent label from the dataset.The second measure metrics is the normalized mutual information (NMI), which is used for determining the quality of clustering.LetC and C
′ denote the sets of cluster from the ground truth and the clustering method, respectively. The mutual information metric M
I
(
C
,
C
′
) is defined as follows:(37)
M
I
C
,
C
′
=
∑
c
i
∈
C
,
c
j
′
∈
C
′
p
c
j
,
c
j
′
·
l
o
g
2
p
c
j
,
c
j
′
p
c
j
·
p
c
j
′
,where p
c
j and p
c
j
′ denote the extent of a sample belonging to the clusters c
i and c
j
′, respectively. p
c
j
,
c
j
′ is a joint probability that an arbitrarily selected sample belongs to c
i as well as c
j
′ at the same time. The NMI can be defined as follows:(38)
N
M
I
C
,
C
′
=
M
I
C
,
C
′
max
⁡
H
C
,
H
C
′
,where H
C and H
C
′ are the entropies of C and C
′, respectively. It can be seen that N
M
I
C
,
C
′ ranges from 0 to 1. The larger the N
M
I is, the better the clustering result will be.


## 4.4. Clustering Results and Analysis

In order to fair examine the clustering performance of all mentioned methods, 20 test runs are conducted on the optimal parameter setting, and then the means and standard deviations are recorded. In each test,Kmeans is applied 20 times with different initialization and the best result is reported. All clustering performances on the three noise contaminated datasets are reported in Tables 3 and 4.Table 3
Clustering performance on local noise contaminated datasets. The best performances are highlighted in bold font.
(a)
Dataset
Accuracy (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
23.10 ± 1.14
26.12 ± 1.27
24.20 ± 1.03
25.05 ± 0.74
27.45 ± 5.22
33.06 ± 1.25
ORL
35.71 ± 2.74
46.99 ± 2.69
44.56 ± 1.59
37.78 ± 1.37
44.78 ± 1.7
54.01 ± 3.3
Yale
32.12 ± 2.38
37.43 ± 2.27
38.18 ± 2.18
32.45 ± 1.16
36.36 ± 1.24
44.42 ± 1.08
(b)
Dataset
Normalized mutual information (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
49.2 ± 2.56
56.85 ± 1.54
54.33 ± 1.12
50.06 ± 2.51
58.26 ± 0.89
67.44 ± 1.73
ORL
53.78 ± 2.56
65.79 ± 1.54
64.57 ± 1.12
58.65 ± 2.51
64.26 ± 0.89
71.50 ± 1.73
Yale
35.73 ± 2.49
42.72 ± 2.30
43.77 ± 1.58
38.84 ± 0.98
42.72 ± 1.12
47.56 ± 1.19Table 4
Clustering performance on global noise contaminated datasets. The best performances are highlighted in bold font.
(a)
Dataset
Accuracy (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
23.76 ± 1.68
27.14 ± 1.57
26.9 ± 1.23
25.76 ± 0.58
28.24 ± 0.73
35.30 ± 0.81
ORL
42.36 ± 3.24
52.78 ± 2.48
52.76 ± 2.85
52.85 ± 0.79
53.91 ± 1.78
61.5 ± 1.92
Yale
36.79 ± 2.24
41.21 ± 2.98
42.03 ± 2.83
40.06 ± 1.43
43.64 ± 1.39
50.30 ± 1.57
(b)
Dataset
Normalized mutual information (%)
Kmeans
PCA
NMF
GNMF
SPNMF
RSPNMF
FERET
50.06 ± 2.74
58.95 ± 1.02
57.94 ± 1.28
60.33 ± 0.60
62.01 ± 0.57
68.65 ± 0.43
ORL
60.14 ± 2.77
70.33 ± 1.39
67.37 ± 1.94
66.7 ± 0.80
70.99 ± 0.74
76.97 ± 0.92
Yale
41.68 ± 2.27
44.78 ± 2.09
46.94 ± 2.37
45.94 ± 0.88
47.68 ± 1.08
54.89 ± 2.08From these results, we can derive the following observations: (a) SPNMF outperformsKmeans, PCA, and NMF in most cases. This demonstrates the importance of structure preservation in discovering the hidden factors. In addition, the performance of SPNMF is superior to GNMF. This reveals that, by leveraging the ratio of local affinity and distant repulsion regularization, SPNMF can learn a better low-dimensional representation. (b) The performance of RSPNMF exceeds SPNMF; this demonstrates the robustness of l
2,1-norm in handling noise data samples.


## 4.5. Parameter Sensitiveness

Our RSPNMF has five essential parameters: the number of nearest neighborsk
1, the number of distant neighbors k
2, the regularization parameter α, the regularization parameter β, and the number of basis vectors r. For all mentioned dimensionality reduction algorithms, the number of basis vectors r is equal to the number of clusters. The complement space and the original space reduced dimension d is set to r
/
2. For all graph-based methods, the number of nearest neighbors is set to {
3,5
,
…
,
19,21
}. For SPNMF and RSPNMF, the number of distant points is also set to {
3,5
,
…
,
19,21
}. For GNMF, the regularization parameter is set by the grid {
10
-
3
,
10
-
2
,
…
,
10
3
}. For SPNMF and RSPNMF, we tune α and β by searching the grid {
10
-
3
,
10
-
2
,
…
,
10
3
}. For SPNMF, γ is set to 1.The result of parameter sensitiveness on different datasets with local and global noise is showed in Figure2. In order to compare with GNMF, we fix the parameter β and tune the parameter α. From Figure 2, we can find that the performance of RSPNMF and SPNMF is superior to GNMF; meanwhile, they are all stable with respect to the change of parameter α. This demonstrates the importance of the distant repulsion preservation on characterizing the intrinsic distribution of data samples. Besides, RSPNMF still outperforms SPNMF when α range within {
0.001,10
}. This is because the robust data reconstruction term of RSPNMF alleviates the noisy data impact on clustering performance.Figure 2
The performance of involved methods with respect to the parameterα on different dataset. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALEThe relevant samples selection plays an important role in all graph regularized NMF method. For GNMF, which only usek
1-nearest neighboring graph to capture the local geometric structure of data distribution; meanwhile, assume that the neighboring samples are with the same semantic label. For noise corrupted data samples, this assumption may not hold. In particular, as k
1 increases, clustering performance of GNMF is more likely to fail as shown in Figure 3. For RSPNMF and SPNMF, the local and distant data samples relationship is also considered. This will definitely reinforce the algorithm robustness to k
1-nearest neighboring samples increase. In addition, from Figure 3 we can find that the introduction of l
2,1-norm improves RSPNMF clustering results.Figure 3
The performance of involved methods with respect to parameterk
1. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALE


## 4.6. Convergence and Computational Cost

The effectiveness of the proposed RSPNMF has been demonstrated by the above clustering experiments. Here, we evaluate the efficiency of RSPNMF, that is, how fast the update rules in (21) and (22) are converged.Figure4 shows the convergence curves of RSPNMF on all the three datasets. For each figure, the y-axis is the value of objective function and the x-axis denotes the iteration number. From the convergence curves, it can be observed that the convergence rate of RSPNMF is fast, and the runtime cost of RSPNMF is acceptable.Figure 4
The convergence curve of RSPNMF with 500 iterations. The result of local noise is shown in the first column. The result of global noise is shown in the second column.
(a)
FERET
(b)
ORL
(c)
YALE


## 5. Conclusion

In this paper, we proposed a novel nonnegative matrix factorization approach named RSPNMF. RSPNMF performs the subspace learning with the nature of manifold structure preservation and robust data reconstruction. As a result, RSPNMF is more discriminating than GNMF which only considers the local structure preservation of the clean data. Experimental results on face image clustering show that RSPNMF provides a better representation in the sense of semantic subspace.The superiority of RSPNMF arises in the following two aspects: (a) the objective function of RSPNMF uses thel
2,1-norm as the discrepancy measure and therefore alleviates the outlier issue that is commonly existing in image clustering. (b) The regularization term of RSPNMF incorporates the local and distant geometric information, which leads the nonnegative matrix factorization to be performed on the specific manifold. The added geometric information generally improves the clustering results.



---

*Source: 7474839-2016-06-29.xml*