# Neighborhood Preserving Convex Nonnegative Matrix Factorization

Authors: Jiang Wei; Li Min; Zhang Yongqing
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The convex nonnegative matrix factorization (CNMF) is a variation of nonnegative matrix factorization (NMF) in which each cluster is expressed by a linear combination of the data points and each data point is represented by a linear combination of the cluster centers. When there exists nonlinearity in the manifold structure, both NMF and CNMF are incapable of characterizing the geometric structure of the data. This paper introduces a neighborhood preserving convex nonnegative matrix factorization (NPCNMF), which imposes an additional constraint on CNMF that each data point can be represented as a linear combination of its neighbors. Thus our method is able to reap the benefits of both nonnegative data factorization and the purpose of manifold structure. An efficient multiplicative updating procedure is produced, and its convergence is guaranteed theoretically. The feasibility and effectiveness of NPCNMF are verified on several standard data sets with promising results.


---

## Body

## 1. Introduction

This nonnegative matrix factorization (NMF) [1, 2] has been widely used in information retrieval, computer vision, pattern recognition, and DNA gene expressions [3, 4]. NMF decomposes the data matrix as the product of two matrices that possess only nonnegative elements. It has been stated by many researchers that there are a lot of favorable properties for such a decomposition over other similar decompositions, such as PCA. One of the most useful properties of NMF is that it usually leads to parts-based representation because it allows only additive, not subtractive, combinations. Such a representation encodes much of the data making them easy to interpret. NMF can be traced back to 1970s and has been studied extensively by Paatero and Tapper [5]. The work of Lee and Seung [1] brought much attention to NMF in machine learning and data mining fields. Since then, various extensions and variations of NMF have been proposed. Li et al. [4] proposed local nonnegative matrix factorization (LNMF) algorithm which imposes extra constraints to the cost function to get more localized and parts-based image features. Hoyer [6, 7] employed sparsity constraints to improve local data representation, while nonnegative tensor factorization was studied in [8, 9] by Hazan et al. to handle the data encoded as high-order tensors. All the methods mentioned above are unsupervised, Wang et al; [10] and Zafeiriou et al. [11] proposed independently the Fisher-NMF, which was further studied by Kotsia et al. [12], by adding an additional constraint seeking to maximize the between-class scatter and minimize the within-class scatter in the subspace spanned by the bases.One of the most important drawbacks of NMF and its variants is the fact that these methods have to be performed in the original feature space of the data points, so that it can not be kernelized and the powerful idea of the kernel method cannot be applied to NMF. Ding et al. [13] proposed convex nonnegative factorization (CNMF) that strives to address the problems while inheriting all the strengths of the above NMF method, which models each cluster as a linear combination of the data points and each data point as a linear combination of the cluster centers. The major advantage of CNMF over NMF is that it can be performed on any data representations, either in the original space or RKHS.Recently, there has been a lot of interest in geometrically motivated approaches to data analysis in high dimensional spaces. When the data lives on or close to a nonlinear low dimensional manifold which is embedded in the high dimensional ambient space [14, 15], Euclidean distance is incapable of charactering the geometric structure of the data and hence traditional methods like NMF and CNMF no longer work well. Both CNMF and NMF do not exploit the geometric structure of the data, which assume that the data points are sampled from a Euclidean space. To address this problem, Cai et al. proposed a graph regularized NMF (GNMF) [16] and locally consistent concept factorization (LCCF) [17], which assume that the nearby data points are likely to be in the same cluster, that is, cluster assumption [18, 19]. The Euclidean and manifold geometry are unified through a regularization framework, which has a better interpretation from manifold perspective.In this paper, we introduce a novel matrix factorization algorithm, called neighborhood preserving convex nonnegative matrix factorization (NPCNMF) which is based on the assumption that if a data point can be reconstructed from its neighbors in the input space, then it can be reconstructed from its neighbors by the same reconstruction coefficients in the low dimensional subspace, that is, local linear embedding assumption [20]. NPCNMF not only inherits the advantages of CNMF, for example, nonnegativity, but also overcomes the shortcomings of CNMF, that is, Euclidean assumption. We also propose a multiplicative algorithm to efficiently solve the corresponding optimization problem and its convergence is theoretically guaranteed.The rest of this paper is organized as follows. In Section2, we briefly review NMF and CNMF. In Section 3, we introduce our algorithm and provide the proof of the convergence of the proposed algorithm. Experiments on three benchmark face recognition data sets are demonstrated in Section 4. Finally, we draw a conclusion and provide suggestions for future work.


## 2. A Review of NMF and CNMF

Nonnegative matrix factorization (NMF) factorizes the data matrix into one nonnegative basis matrix and one nonnegative coefficient matrix. Given a nonnegative dataX
=
[
x
1
,
x
2
,
…
,
x
N
]
∈
ℝ
+
m
×
N, each column of X is a sample point. NMF aims to find two nonnegative matrices U
∈
ℝ
+
m
×
r and V
∈
ℝ
+
N
×
r which minimize the following objective function:
(1)
𝒥
NMF
=
∥
X
-
U
V
T
∥
2
s
.
t
.
U
≥
0
,
V
≥
0
,
where ∥
·
∥
F is Frobenius norm.The objective function is joint optimization problem of basis matrixU and coefficient matrix V. Although it is not jointly convex to U and V, it is convex with respect to each of them when the other one is fixed. Therefore, it is unrealistic to expect an algorithm to find the global minimum of 𝒥
NMF. To optimize the objective, Lee and Seung [2] presented an iterative multiplicative updating algorithm as follows:
(2)
u
i
k
t
+
1
=
u
i
k
t
(
X
V
)
i
k
(
U
V
T
V
)
i
k
,
v
j
k
t
+
1
=
v
j
k
t
(
X
T
U
)
j
k
(
V
U
T
U
)
j
k
.It is proved that the above updated steps will find a local minimum of the objective function in (1).In reality, we haver
≪
m and r
≪
N. Thus, NMF essentially tries to find a compressed approximation of the original matrix, X
≈
U
V
T. We can view this approximation column as follows:
(3)
x
j
=
∑
k
=
1
r
‍
u
k
v
j
k
,
where u
k is the kth column vector of U. Thus, each data vector x
j is approximated by a linear combination of the columns of U, weighted by the components of V. One limitation of NMF is that the nonnegative requirement is not applicable to applications where the data involves negative number. The second is that it is not clear how to effectively perform NMF in the transformed data space so that the powerful kernel method can be applied. To overcome the problem, Ding et al. [13] proposed a convex nonnegative matrix factorization (CNMF) algorithm where nonnegative and mixed-sign data matrices are applied. CNMF claims that each base can be characterized by a linear combination of the entire data points while each data point can be approximated by a linear combination of all the bases. Translating the statements into mathematics, we have
(4)
u
k
=
∑
j
=
1
N
‍
x
j
w
j
k
,
(5)
x
j
=
∑
k
=
1
r
‍
u
k
v
j
k
,
where w
j
k is a nonnegative weight in which data point x
j is related to kth base and v
j
k is a nonnegative projection value of x
j. Replacing u
k in (5) with (4), we have
(6)
x
j
=
∑
k
=
1
r
‍
∑
j
=
1
N
‍
x
j
w
j
k
u
j
k
.
We form the m
×
N data matrix X
=
[
x
1
,
x
2
,
…
,
x
N
] using the feature vector of data point x
i as the ith column, the N
×
r matrix W
=
[
w
i
j
] using bases w
i
j, and N
×
r projection matrix V
=
[
v
i
j
] using the projection values v
i
j. From (6), we have
(7)
X
=
X
W
V
T
.Equation (7) can be interpreted as the approximation of the original data set. Minimizing the squared error 𝒥 and its approximation [13]
(8)
𝒥
CNMF
=
∥
X
-
X
W
V
T
∥
2
,
where X
∈
ℝ
m
×
N, W
∈
ℝ
N
×
r, V
∈
ℝ
N
×
r. The matrices W and V are updated iteratively until convergence using the following rules:
(9)
w
i
k
=
w
i
k
(
Y
+
V
)
i
k
+
(
Y
-
W
V
T
V
)
i
k
(
Y
-
V
)
i
k
+
(
Y
+
W
V
T
V
)
i
k
,
v
i
k
=
v
i
k
(
Y
+
W
)
i
k
+
(
V
W
T
Y
-
W
)
i
k
(
Y
-
W
)
i
k
+
(
V
W
T
Y
+
W
)
i
k
,
where Y
=
X
T
X and the matrix Y
+ and Y
- are given by
(10)
Y
i
k
+
=
1
2
(
|
Y
i
k
|
+
Y
i
k
)
,
Y
i
k
-
=
1
2
(
|
Y
i
k
|
-
Y
i
k
)
,
respectively.


## 3. Neighborhood Preserving Convex Nonnegative Matrix Factorization

In this section, we introduce our neighborhood preserving convex nonnegative matrix factorization method, which takes the local linear embedding constraint as an additional requirement. The method presented in this paper is fundamentally motivated from the neighborhood preserving embedding.


### 3.1. The Objective Function

Many real world data are actually sampled from a nonlinear low dimensional manifold which is embedded in the high dimensional ambient space. Both NMF and CF perform the factorization in the Euclidean space. They fail to discover the local geometrical structure of the data space, which is essential to the clustering problem. NPE aims at preserving the local manifold structure. Specifically, for each data point, it is represented as a linear combination of the neighboring data points and the combination coefficients are specified in the weight matrix. We can find an optimal embedding such that the combination coefficients can be preserved in the low dimensional subspace.For each data point, we find itsk nearest neighbors. And we can characterize the local geometrical structure by linear coefficients that reconstruct each data point from its neighbors. The reconstruction coefficients are computed by the following objective function:
(11)
min
⁡
∥
x
i
-
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
x
j
∥
2
s
.
t
.
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
=
1
,
and w
i
j
=
0 if x
j
∉
𝒩
k
(
x
i
), where 𝒩
k
(
x
i
) denotes the k nearest neighborhood of x
i.Thenv
i, 1
≤
i
≤
n in the dimensionality reduced space can be preserved by minimizing
(12)
min
⁡
⁡
∥
v
i
-
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
v
j
∥
2
=
tr
⁡
(
V
T
(
I
-
W
)
(
I
-
W
)
V
)
=
tr
⁡
(
V
T
L
V
)
,
where tr
⁡
(
·
) denotes the trace of a matrix, I
∈
ℝ
n
×
n is an identity matrix, and L
=
(
I
-
W
)
(
I
-
W
). By minimizing (12), we essentially try to formalize our intuition that if a data point can be represented from its neighbors in the original space, then it can be represented from its neighbors by the same combination coefficients in the dimensionality reduced space.With the neighborhood preserving constraint, CNMF incorporates (8) and minimizes the objective function as follows:
(13)
𝒥
=
∥
X
-
X
W
V
T
∥
2
+
λ
tr
⁡
(
V
T
L
V
)
,
where λ is a positive regularization parameter controlling the contribution of the additional constraint. We call (13) neighborhood preserving convex nonnegative matrix factorization (NPCNMF). Let λ
=
0; (13) degenerates to the original CNMF.


### 3.2. The Algorithm

We introduce an iterative algorithm to find a local minimum for the optimization problem. By definingK
=
X
T
X and using the matrix properties ∥
A
∥
2
=
tr
⁡
(
A
T
A
), tr
⁡
(
A
B
)
=
tr
⁡
(
B
A
), and tr
⁡
(
A
)
=
tr
⁡
(
A
T
), we can rewrite the objective function 𝒥 as follows:
(14)
𝒥
=
tr
⁡
(
(
X
-
X
W
V
T
)
T
(
X
-
X
W
V
T
)
)
+
λ
tr
⁡
(
V
T
L
V
)
=
tr
⁡
(
(
I
-
W
V
T
)
T
K
(
I
-
W
V
T
)
)
+
λ
tr
⁡
(
V
T
L
V
)
=
tr
⁡
(
K
)
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
.This is a typical constrained optimization problem and can be solved using the Lagrange multiplier method. Letψ
i
j and ϕ
i
j be the Lagrange multiplier for constraint w
i
j
≥
0 and ν
i
j
≥
0, respectively, and let Ψ
=
[
ψ
i
j
] and Φ
=
[
ϕ
i
j
]. The Lagrangian function is
(15)
ℒ
=
tr
⁡
(
K
)
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
+
tr
⁡
(
Ψ
W
T
)
+
tr
⁡
(
Φ
V
T
)
.The partial derivatives ofℒ with respect to W and V are
(16)
∂
ℒ
∂
W
=
-
2
K
V
+
2
K
W
V
T
V
+
Ψ
,
∂
ℒ
∂
V
=
-
2
K
W
+
2
V
W
T
K
W
+
2
λ
L
V
+
Φ
.Using the Karush-Kuhn-Tucker conditionsψ
i
k
w
i
k
=
0 and ϕ
j
k
ν
j
k
=
0, we get the following equations for w
i
k and v
j
k:
(17)
(
-
K
V
+
K
W
V
T
V
)
i
k
w
i
k
=
0
,
(
-
K
W
+
V
W
T
K
W
+
λ
L
V
)
j
k
v
j
k
=
0
.The corresponding equivalent formulas are as follows:(18)
(
-
K
V
+
K
W
V
T
V
)
i
k
w
i
k
2
=
0
,
(
-
K
W
+
V
W
T
K
W
+
λ
L
V
)
j
k
v
j
k
2
=
0
.Introduce(19)
A
=
A
+
-
A
-
,
where A
i
k
+
=
(
|
A
i
k
|
+
A
i
k
)
/
2 and A
i
k
-
=
(
|
A
i
k
|
-
A
i
k
)
/
2.The equations lead to the following updating formulas:(20)
w
i
k
=
w
i
k
(
K
+
V
)
i
k
+
(
K
-
W
V
T
V
)
i
k
(
K
-
V
)
i
k
+
(
K
+
W
V
T
V
)
i
k
,
(21)
v
j
k
=
v
i
k
(
K
+
W
)
j
k
+
(
V
W
T
K
-
W
)
j
k
+
λ
(
L
-
V
)
j
k
(
K
-
W
)
j
k
+
(
V
W
T
K
+
W
)
j
k
+
λ
(
L
+
V
)
j
k
.Note that the solution to minimizing the criterion function𝒥 is not unique. If W and V are the solution, then, W
D, W
D
-
1 will also form a solution for any positive diagonal matrix D. To make the solution unique, we will further require that w
T
K
w
=
1, where w is the column vector of W. The matrix V will be adjusted accordingly so that W
V
T does not change. This can be achieved by
(22)
V
⟵
V
[
diag
⁡
(
W
T
K
W
)
]
1
/
2
,
W
⟵
W
[
diag
⁡
(
W
T
K
W
)
]
1
/
2
.


### 3.3. Convergence Analysis

In this section, we will investigate the convergence of the updating formula in (14). We use the auxiliary function approach [16] to prove the convergence. Here we first introduce the definition of auxiliary function [16].Definition 1.
Z
(
h
,
h
′
) is an auxiliary function of F
(
h
) if the conditions
(23)
Z
(
h
,
h
′
)
≥
F
(
h
)
,
Z
(
h
,
h
)
=
F
(
h
)
are satisfied.Lemma 2.
IfZ is an auxiliary function for F, then F is nonincreasing under the update
(24)
h
(
t
+
1
)
=
arg
min
⁡
h
Z
(
h
,
h
(
t
)
)
.Proof.
Consider(25)
F
(
h
(
t
+
1
)
)
≤
Z
(
h
(
t
+
1
)
,
h
(
t
)
)
≤
Z
(
h
(
t
)
,
h
(
t
)
)
=
F
(
h
(
t
)
)
.Lemma 3.
For any nonnegative matricesA
∈
ℝ
n
×
n, B
∈
ℝ
k
×
k, S
∈
ℝ
n
×
k, and  S
′
∈
ℝ
n
×
k, and A, B are symmetric, then the following inequality holds:
(26)
∑
i
=
1
n
‍
∑
p
=
1
k
‍
(
A
S
′
B
)
i
p
2
S
i
p
′
≥
tr
⁡
(
S
T
A
S
B
)
.The correctness and convergence of the algorithm are addressed in the following.For givenK, fixing V, considering any element w
i
j in W, we use 𝒥
(
W
) to denote the part of 𝒥, which is only relevant to w
i
j. We get
(27)
𝒥
(
W
)
=
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
.Theorem 4.
One rewrites𝒥
(
W
) as follows:
(28)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
T
A
+
H
C
a
m
r
o
-
H
T
A
-
H
C
)
,
where B
=
K
V, A
=
K, C
=
V
T
V, and H
=
W.
Then the following function(29)
Z
(
H
,
H
′
)
=
-
∑
i
k
‍
2
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
+
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
H
i
k
′
+
∑
i
k
‍
(
A
+
H
′
C
)
i
k
H
i
k
2
H
i
k
′
-
∑
i
j
k
l
‍
A
i
j
-
H
j
k
′
C
k
l
H
i
l
′
(
1
+
log
⁡
H
j
k
H
i
l
H
j
k
′
H
i
l
′
)
is an auxiliary function of 𝒥
(
H
): that is, it satisfies the requirements 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
). Furthermore, it is a convex function of H and its global minimum is
(30)
H
i
k
=
argmin
H
H
i
k
′
B
i
k
+
+
(
H
′
A
-
C
)
i
k
B
i
k
-
+
(
A
+
H
′
C
)
i
k
.
From its minima and settingH
(
t
+
1
)
←
H and H
(
t
)
←
H
′, one recovers (20), letting B
+
=
(
K
)
+
V, B
-
=
(
K
)
-
V, A
=
K, C
=
V
T
V, and H
=
W.Proof.
The function𝒥
(
H
) is
(31)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
T
A
+
H
C
-
H
T
A
-
H
C
)
.
We find upper bounds for each of the two positive terms and lower bounds for each of the two negative terms. For the third term in𝒥
(
H
), by applying Lemma 3, we obtain an upper bound
(32)
tr
⁡
(
H
T
A
+
H
C
)
≤
∑
i
k
‍
(
A
+
H
′
C
)
i
k
H
i
k
2
H
i
k
′
.
The second term of𝒥
(
H
) is bounded by
(33)
tr
⁡
(
H
T
B
-
)
=
∑
i
k
‍
H
i
k
B
i
k
-
≤
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
2
H
i
k
′
.
To obtain lower bounds for the two remaining terms, we use the inequalityz
≥
1
+
log
⁡
z, which holds for any z
>
0, and the first term in 𝒥
(
H
) is bounded by
(34)
tr
⁡
(
H
T
B
+
)
=
∑
i
k
‍
B
i
k
+
H
i
k
≥
∑
i
k
‍
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
.
The last term in𝒥
(
H
) is bounded by
(35)
tr
⁡
(
H
T
A
-
H
C
)
≥
∑
i
j
k
l
‍
A
i
j
-
H
j
k
′
C
k
l
H
i
l
′
(
1
+
log
⁡
H
j
k
H
i
l
H
j
k
′
H
i
l
′
)
.
Collecting all bounds, we obtainZ
(
H
,
H
′
) as in (29). Obviously, 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
).
To find the minimum ofZ
(
H
,
H
′
), we take
(36)
∂
Z
(
H
,
H
′
)
∂
H
i
k
=
-
2
B
i
k
+
H
i
k
′
H
i
k
+
2
B
i
k
-
H
i
k
H
i
k
′
+
2
(
H
′
A
+
C
)
i
k
H
i
k
H
i
k
′
-
2
(
H
′
A
-
C
)
i
k
H
i
k
′
H
i
k
.
To find the minimum ofZ
(
H
,
H
′
), we take the Hessian matrix of Z
(
H
,
H
′
)
(37)
∂
2
Z
(
H
,
H
′
)
∂
H
i
k
∂
H
j
l
=
δ
i
j
δ
k
l
Y
i
k
To be a diagonal matrix with positive entries
(38)
Y
i
k
=
4
[
(
B
+
)
i
k
+
(
H
′
A
-
C
)
i
k
]
H
i
k
′
H
i
k
2
+
2
B
i
k
-
+
(
H
′
A
+
C
)
i
k
H
i
k
′
.
Thus,Z
(
H
,
H
′
) is a convex function of H. Therefore, we obtain the global minimum by setting ∂
Z
(
H
,
H
′
)
/
∂
H
i
k
=
0 in (36) and solving for H. Rearranging, we obtain (30).Theorem 5.
UpdatingW using (20) will monotonically decrease the value of the objective in (13); hence it converges.Proof.
By Lemma2 and Theorem 4, we can get that 𝒥
(
W
t
)
=
Z
(
W
(
t
)
,
W
(
t
)
)
≥
Z
(
W
(
t
+
1
)
,
W
(
t
)
)
≥
𝒥
(
W
t
+
1
), so 𝒥
(
W
) is monotonically decreasing. Since 𝒥
(
W
) is obviously bounded below, we prove this theorem.
For givenK, fixing W, considering any element v
i
j in V, we use 𝒥
(
V
) to denote the part of 𝒥, which is only relevant to v
i
j. We get
(39)
𝒥
(
V
)
=
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
.Theorem 6.
One rewrites𝒥
(
V
) as follows:
(40)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
A
+
H
T
-
H
A
-
H
T
+
λ
H
T
L
+
H
-
λ
H
T
L
-
H
)
,
where B
=
K
W, A
=
W
T
X
T
X
W, and H
=
V.
Then the following function(41)
Z
(
H
,
H
′
)
=
-
∑
i
k
‍
2
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
+
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
H
i
k
′
+
∑
i
k
‍
(
A
+
H
′
)
i
k
H
i
k
2
H
i
k
′
-
∑
i
k
l
‍
A
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
+
λ
∑
i
k
‍
(
L
+
H
′
)
i
k
H
i
k
2
H
i
k
′
-
λ
∑
i
k
l
‍
L
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
is an auxiliary function of 𝒥
(
H
): that is, it satisfies the requirements 𝒥
(
H
)
≤
Z
(
H
,
H
) and 𝒥
(
H
)
=
Z
(
H
,
H
′
). Furthermore, it is a convex function of H and its global minimum is
(42)
H
i
k
=
H
i
k
′
B
i
k
+
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
.
From its minima and settingH
(
t
+
1
)
←
H and H
(
t
)
←
H
′, one recovers (21), letting B
+
=
(
K
)
+
U, B
-
=
(
K
)
-
U, A
=
W
T
K
W, and H
=
V.Proof.
The function𝒥
(
H
) is
(43)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
A
+
H
T
-
H
A
-
H
T
+
λ
H
T
L
+
H
-
λ
H
T
L
-
H
)
.
We find upper bounds for each of the three positive terms and lower bounds for each of the three negative terms. For the third term in𝒥
(
H
), by applying Lemma 3 and setting A
←
I, B
←
A
+, we obtain an upper bound
(44)
tr
⁡
(
H
A
+
H
T
)
≤
∑
i
k
‍
(
H
′
A
+
)
i
k
H
i
k
2
H
i
k
′
.
The second term of𝒥
(
H
) is bounded by
(45)
tr
⁡
(
H
T
B
-
)
=
∑
i
k
‍
H
i
k
B
i
k
-
≤
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
2
H
i
k
′
,
using the inequality a
≤
(
a
2
+
b
2
)
/
2
b, which holds for any a
,
b
>
0.
For the fifth term in𝒥
(
H
), setting A
←
L
+, B
←
I, and S
←
H, we obtain an upper bound
(46)
λ
tr
⁡
(
H
T
L
+
H
)
≤
λ
∑
i
k
‍
(
L
+
H
′
)
i
k
H
i
k
2
H
i
k
′
.
To obtain lower bounds for the three remaining terms, we use the inequalityz
≥
1
+
log
⁡
z, which holds for any z
>
0, and the first term in J
(
H
) is bounded by
(47)
tr
⁡
(
H
T
B
+
)
=
∑
i
k
‍
B
i
k
+
H
i
k
≥
∑
i
k
‍
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
.
The fourth term in𝒥
(
H
) is bounded by
(48)
tr
⁡
(
H
A
-
H
T
)
≥
∑
i
k
l
‍
A
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
.
The last term in𝒥
(
H
) is bounded by
(49)
λ
tr
⁡
(
H
T
L
-
H
)
≥
λ
∑
i
k
l
‍
L
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
.
Collecting all bounds, we obtainZ
(
H
,
H
′
) as in (41). Obviously, 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
).
To find the minimum ofZ
(
H
,
H
′
), we take
(50)
∂
Z
(
H
,
H
′
)
∂
H
i
k
=
-
2
B
i
k
+
H
i
k
′
H
i
k
+
2
B
i
k
-
H
i
k
H
i
k
′
+
2
(
H
′
A
+
)
i
k
H
i
k
H
i
k
′
-
2
(
H
′
A
-
)
i
k
H
i
k
′
H
i
k
+
2
λ
(
L
+
H
′
)
i
k
H
i
k
H
i
k
′
-
2
λ
(
L
-
H
′
)
i
k
H
i
k
′
H
i
k
=
0
.
We have(51)
-
2
B
i
k
+
(
H
′
)
i
k
2
+
2
B
i
k
-
H
i
k
2
+
2
(
H
′
A
+
)
i
k
H
i
k
2
-
2
(
H
′
A
-
)
i
k
(
H
′
)
i
k
2
+
2
λ
(
L
+
H
′
)
i
k
H
i
k
2
-
2
λ
(
L
-
H
′
)
i
k
(
H
′
)
i
k
2
=
0
.
Therefore(52)
H
i
k
=
H
i
k
′
B
i
k
+
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
.
The Hessian matrix containing the second derivatives(53)
∂
2
Z
(
H
,
H
′
)
∂
H
i
k
∂
H
j
l
=
δ
i
j
δ
k
l
Y
i
k
is a diagonal matrix with positive entries
(54)
Y
i
k
=
4
[
(
B
+
)
i
k
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
]
H
i
k
′
H
i
k
2
+
2
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
H
i
k
′
.
Thus,Z
(
H
,
H
′
) is a convex function of H. Therefore, we obtain the global minimum by setting ∂
Z
(
H
,
H
′
)
/
∂
H
i
k
=
0 in (41) and solving for H. Rearranging, we obtain  (21).Theorem 7.
UpdatingV using (21) will monotonically decrease the value of the objective in (13); hence it converges.Proof.
By Lemma2 and Theorem 4, we can get that 𝒥
(
V
t
)
=
Z
(
V
(
t
)
,
V
(
t
)
)
≥
Z
(
V
(
t
+
1
)
,
V
(
t
)
)
≥
𝒥
(
V
t
+
1
), so 𝒥
(
V
) is monotonically decreasing. Since 𝒥
(
V
) is obviously bounded below, we prove this theorem.


## 3.1. The Objective Function

Many real world data are actually sampled from a nonlinear low dimensional manifold which is embedded in the high dimensional ambient space. Both NMF and CF perform the factorization in the Euclidean space. They fail to discover the local geometrical structure of the data space, which is essential to the clustering problem. NPE aims at preserving the local manifold structure. Specifically, for each data point, it is represented as a linear combination of the neighboring data points and the combination coefficients are specified in the weight matrix. We can find an optimal embedding such that the combination coefficients can be preserved in the low dimensional subspace.For each data point, we find itsk nearest neighbors. And we can characterize the local geometrical structure by linear coefficients that reconstruct each data point from its neighbors. The reconstruction coefficients are computed by the following objective function:
(11)
min
⁡
∥
x
i
-
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
x
j
∥
2
s
.
t
.
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
=
1
,
and w
i
j
=
0 if x
j
∉
𝒩
k
(
x
i
), where 𝒩
k
(
x
i
) denotes the k nearest neighborhood of x
i.Thenv
i, 1
≤
i
≤
n in the dimensionality reduced space can be preserved by minimizing
(12)
min
⁡
⁡
∥
v
i
-
∑
x
j
∈
𝒩
k
(
x
i
)
‍
w
i
j
v
j
∥
2
=
tr
⁡
(
V
T
(
I
-
W
)
(
I
-
W
)
V
)
=
tr
⁡
(
V
T
L
V
)
,
where tr
⁡
(
·
) denotes the trace of a matrix, I
∈
ℝ
n
×
n is an identity matrix, and L
=
(
I
-
W
)
(
I
-
W
). By minimizing (12), we essentially try to formalize our intuition that if a data point can be represented from its neighbors in the original space, then it can be represented from its neighbors by the same combination coefficients in the dimensionality reduced space.With the neighborhood preserving constraint, CNMF incorporates (8) and minimizes the objective function as follows:
(13)
𝒥
=
∥
X
-
X
W
V
T
∥
2
+
λ
tr
⁡
(
V
T
L
V
)
,
where λ is a positive regularization parameter controlling the contribution of the additional constraint. We call (13) neighborhood preserving convex nonnegative matrix factorization (NPCNMF). Let λ
=
0; (13) degenerates to the original CNMF.


## 3.2. The Algorithm

We introduce an iterative algorithm to find a local minimum for the optimization problem. By definingK
=
X
T
X and using the matrix properties ∥
A
∥
2
=
tr
⁡
(
A
T
A
), tr
⁡
(
A
B
)
=
tr
⁡
(
B
A
), and tr
⁡
(
A
)
=
tr
⁡
(
A
T
), we can rewrite the objective function 𝒥 as follows:
(14)
𝒥
=
tr
⁡
(
(
X
-
X
W
V
T
)
T
(
X
-
X
W
V
T
)
)
+
λ
tr
⁡
(
V
T
L
V
)
=
tr
⁡
(
(
I
-
W
V
T
)
T
K
(
I
-
W
V
T
)
)
+
λ
tr
⁡
(
V
T
L
V
)
=
tr
⁡
(
K
)
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
.This is a typical constrained optimization problem and can be solved using the Lagrange multiplier method. Letψ
i
j and ϕ
i
j be the Lagrange multiplier for constraint w
i
j
≥
0 and ν
i
j
≥
0, respectively, and let Ψ
=
[
ψ
i
j
] and Φ
=
[
ϕ
i
j
]. The Lagrangian function is
(15)
ℒ
=
tr
⁡
(
K
)
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
+
tr
⁡
(
Ψ
W
T
)
+
tr
⁡
(
Φ
V
T
)
.The partial derivatives ofℒ with respect to W and V are
(16)
∂
ℒ
∂
W
=
-
2
K
V
+
2
K
W
V
T
V
+
Ψ
,
∂
ℒ
∂
V
=
-
2
K
W
+
2
V
W
T
K
W
+
2
λ
L
V
+
Φ
.Using the Karush-Kuhn-Tucker conditionsψ
i
k
w
i
k
=
0 and ϕ
j
k
ν
j
k
=
0, we get the following equations for w
i
k and v
j
k:
(17)
(
-
K
V
+
K
W
V
T
V
)
i
k
w
i
k
=
0
,
(
-
K
W
+
V
W
T
K
W
+
λ
L
V
)
j
k
v
j
k
=
0
.The corresponding equivalent formulas are as follows:(18)
(
-
K
V
+
K
W
V
T
V
)
i
k
w
i
k
2
=
0
,
(
-
K
W
+
V
W
T
K
W
+
λ
L
V
)
j
k
v
j
k
2
=
0
.Introduce(19)
A
=
A
+
-
A
-
,
where A
i
k
+
=
(
|
A
i
k
|
+
A
i
k
)
/
2 and A
i
k
-
=
(
|
A
i
k
|
-
A
i
k
)
/
2.The equations lead to the following updating formulas:(20)
w
i
k
=
w
i
k
(
K
+
V
)
i
k
+
(
K
-
W
V
T
V
)
i
k
(
K
-
V
)
i
k
+
(
K
+
W
V
T
V
)
i
k
,
(21)
v
j
k
=
v
i
k
(
K
+
W
)
j
k
+
(
V
W
T
K
-
W
)
j
k
+
λ
(
L
-
V
)
j
k
(
K
-
W
)
j
k
+
(
V
W
T
K
+
W
)
j
k
+
λ
(
L
+
V
)
j
k
.Note that the solution to minimizing the criterion function𝒥 is not unique. If W and V are the solution, then, W
D, W
D
-
1 will also form a solution for any positive diagonal matrix D. To make the solution unique, we will further require that w
T
K
w
=
1, where w is the column vector of W. The matrix V will be adjusted accordingly so that W
V
T does not change. This can be achieved by
(22)
V
⟵
V
[
diag
⁡
(
W
T
K
W
)
]
1
/
2
,
W
⟵
W
[
diag
⁡
(
W
T
K
W
)
]
1
/
2
.


## 3.3. Convergence Analysis

In this section, we will investigate the convergence of the updating formula in (14). We use the auxiliary function approach [16] to prove the convergence. Here we first introduce the definition of auxiliary function [16].Definition 1.
Z
(
h
,
h
′
) is an auxiliary function of F
(
h
) if the conditions
(23)
Z
(
h
,
h
′
)
≥
F
(
h
)
,
Z
(
h
,
h
)
=
F
(
h
)
are satisfied.Lemma 2.
IfZ is an auxiliary function for F, then F is nonincreasing under the update
(24)
h
(
t
+
1
)
=
arg
min
⁡
h
Z
(
h
,
h
(
t
)
)
.Proof.
Consider(25)
F
(
h
(
t
+
1
)
)
≤
Z
(
h
(
t
+
1
)
,
h
(
t
)
)
≤
Z
(
h
(
t
)
,
h
(
t
)
)
=
F
(
h
(
t
)
)
.Lemma 3.
For any nonnegative matricesA
∈
ℝ
n
×
n, B
∈
ℝ
k
×
k, S
∈
ℝ
n
×
k, and  S
′
∈
ℝ
n
×
k, and A, B are symmetric, then the following inequality holds:
(26)
∑
i
=
1
n
‍
∑
p
=
1
k
‍
(
A
S
′
B
)
i
p
2
S
i
p
′
≥
tr
⁡
(
S
T
A
S
B
)
.The correctness and convergence of the algorithm are addressed in the following.For givenK, fixing V, considering any element w
i
j in W, we use 𝒥
(
W
) to denote the part of 𝒥, which is only relevant to w
i
j. We get
(27)
𝒥
(
W
)
=
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
.Theorem 4.
One rewrites𝒥
(
W
) as follows:
(28)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
T
A
+
H
C
a
m
r
o
-
H
T
A
-
H
C
)
,
where B
=
K
V, A
=
K, C
=
V
T
V, and H
=
W.
Then the following function(29)
Z
(
H
,
H
′
)
=
-
∑
i
k
‍
2
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
+
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
H
i
k
′
+
∑
i
k
‍
(
A
+
H
′
C
)
i
k
H
i
k
2
H
i
k
′
-
∑
i
j
k
l
‍
A
i
j
-
H
j
k
′
C
k
l
H
i
l
′
(
1
+
log
⁡
H
j
k
H
i
l
H
j
k
′
H
i
l
′
)
is an auxiliary function of 𝒥
(
H
): that is, it satisfies the requirements 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
). Furthermore, it is a convex function of H and its global minimum is
(30)
H
i
k
=
argmin
H
H
i
k
′
B
i
k
+
+
(
H
′
A
-
C
)
i
k
B
i
k
-
+
(
A
+
H
′
C
)
i
k
.
From its minima and settingH
(
t
+
1
)
←
H and H
(
t
)
←
H
′, one recovers (20), letting B
+
=
(
K
)
+
V, B
-
=
(
K
)
-
V, A
=
K, C
=
V
T
V, and H
=
W.Proof.
The function𝒥
(
H
) is
(31)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
T
A
+
H
C
-
H
T
A
-
H
C
)
.
We find upper bounds for each of the two positive terms and lower bounds for each of the two negative terms. For the third term in𝒥
(
H
), by applying Lemma 3, we obtain an upper bound
(32)
tr
⁡
(
H
T
A
+
H
C
)
≤
∑
i
k
‍
(
A
+
H
′
C
)
i
k
H
i
k
2
H
i
k
′
.
The second term of𝒥
(
H
) is bounded by
(33)
tr
⁡
(
H
T
B
-
)
=
∑
i
k
‍
H
i
k
B
i
k
-
≤
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
2
H
i
k
′
.
To obtain lower bounds for the two remaining terms, we use the inequalityz
≥
1
+
log
⁡
z, which holds for any z
>
0, and the first term in 𝒥
(
H
) is bounded by
(34)
tr
⁡
(
H
T
B
+
)
=
∑
i
k
‍
B
i
k
+
H
i
k
≥
∑
i
k
‍
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
.
The last term in𝒥
(
H
) is bounded by
(35)
tr
⁡
(
H
T
A
-
H
C
)
≥
∑
i
j
k
l
‍
A
i
j
-
H
j
k
′
C
k
l
H
i
l
′
(
1
+
log
⁡
H
j
k
H
i
l
H
j
k
′
H
i
l
′
)
.
Collecting all bounds, we obtainZ
(
H
,
H
′
) as in (29). Obviously, 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
).
To find the minimum ofZ
(
H
,
H
′
), we take
(36)
∂
Z
(
H
,
H
′
)
∂
H
i
k
=
-
2
B
i
k
+
H
i
k
′
H
i
k
+
2
B
i
k
-
H
i
k
H
i
k
′
+
2
(
H
′
A
+
C
)
i
k
H
i
k
H
i
k
′
-
2
(
H
′
A
-
C
)
i
k
H
i
k
′
H
i
k
.
To find the minimum ofZ
(
H
,
H
′
), we take the Hessian matrix of Z
(
H
,
H
′
)
(37)
∂
2
Z
(
H
,
H
′
)
∂
H
i
k
∂
H
j
l
=
δ
i
j
δ
k
l
Y
i
k
To be a diagonal matrix with positive entries
(38)
Y
i
k
=
4
[
(
B
+
)
i
k
+
(
H
′
A
-
C
)
i
k
]
H
i
k
′
H
i
k
2
+
2
B
i
k
-
+
(
H
′
A
+
C
)
i
k
H
i
k
′
.
Thus,Z
(
H
,
H
′
) is a convex function of H. Therefore, we obtain the global minimum by setting ∂
Z
(
H
,
H
′
)
/
∂
H
i
k
=
0 in (36) and solving for H. Rearranging, we obtain (30).Theorem 5.
UpdatingW using (20) will monotonically decrease the value of the objective in (13); hence it converges.Proof.
By Lemma2 and Theorem 4, we can get that 𝒥
(
W
t
)
=
Z
(
W
(
t
)
,
W
(
t
)
)
≥
Z
(
W
(
t
+
1
)
,
W
(
t
)
)
≥
𝒥
(
W
t
+
1
), so 𝒥
(
W
) is monotonically decreasing. Since 𝒥
(
W
) is obviously bounded below, we prove this theorem.
For givenK, fixing W, considering any element v
i
j in V, we use 𝒥
(
V
) to denote the part of 𝒥, which is only relevant to v
i
j. We get
(39)
𝒥
(
V
)
=
-
2
tr
⁡
(
V
W
T
K
)
+
tr
⁡
(
V
W
T
K
W
V
T
)
+
λ
tr
⁡
(
V
T
L
V
)
.Theorem 6.
One rewrites𝒥
(
V
) as follows:
(40)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
A
+
H
T
-
H
A
-
H
T
+
λ
H
T
L
+
H
-
λ
H
T
L
-
H
)
,
where B
=
K
W, A
=
W
T
X
T
X
W, and H
=
V.
Then the following function(41)
Z
(
H
,
H
′
)
=
-
∑
i
k
‍
2
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
+
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
H
i
k
′
+
∑
i
k
‍
(
A
+
H
′
)
i
k
H
i
k
2
H
i
k
′
-
∑
i
k
l
‍
A
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
+
λ
∑
i
k
‍
(
L
+
H
′
)
i
k
H
i
k
2
H
i
k
′
-
λ
∑
i
k
l
‍
L
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
is an auxiliary function of 𝒥
(
H
): that is, it satisfies the requirements 𝒥
(
H
)
≤
Z
(
H
,
H
) and 𝒥
(
H
)
=
Z
(
H
,
H
′
). Furthermore, it is a convex function of H and its global minimum is
(42)
H
i
k
=
H
i
k
′
B
i
k
+
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
.
From its minima and settingH
(
t
+
1
)
←
H and H
(
t
)
←
H
′, one recovers (21), letting B
+
=
(
K
)
+
U, B
-
=
(
K
)
-
U, A
=
W
T
K
W, and H
=
V.Proof.
The function𝒥
(
H
) is
(43)
𝒥
(
H
)
=
tr
⁡
(
-
2
H
T
B
+
+
2
H
T
B
-
+
H
A
+
H
T
-
H
A
-
H
T
+
λ
H
T
L
+
H
-
λ
H
T
L
-
H
)
.
We find upper bounds for each of the three positive terms and lower bounds for each of the three negative terms. For the third term in𝒥
(
H
), by applying Lemma 3 and setting A
←
I, B
←
A
+, we obtain an upper bound
(44)
tr
⁡
(
H
A
+
H
T
)
≤
∑
i
k
‍
(
H
′
A
+
)
i
k
H
i
k
2
H
i
k
′
.
The second term of𝒥
(
H
) is bounded by
(45)
tr
⁡
(
H
T
B
-
)
=
∑
i
k
‍
H
i
k
B
i
k
-
≤
∑
i
k
‍
B
i
k
-
H
i
k
2
+
(
H
′
)
i
k
2
2
H
i
k
′
,
using the inequality a
≤
(
a
2
+
b
2
)
/
2
b, which holds for any a
,
b
>
0.
For the fifth term in𝒥
(
H
), setting A
←
L
+, B
←
I, and S
←
H, we obtain an upper bound
(46)
λ
tr
⁡
(
H
T
L
+
H
)
≤
λ
∑
i
k
‍
(
L
+
H
′
)
i
k
H
i
k
2
H
i
k
′
.
To obtain lower bounds for the three remaining terms, we use the inequalityz
≥
1
+
log
⁡
z, which holds for any z
>
0, and the first term in J
(
H
) is bounded by
(47)
tr
⁡
(
H
T
B
+
)
=
∑
i
k
‍
B
i
k
+
H
i
k
≥
∑
i
k
‍
B
i
k
+
H
i
k
′
(
1
+
log
⁡
H
i
k
H
i
k
′
)
.
The fourth term in𝒥
(
H
) is bounded by
(48)
tr
⁡
(
H
A
-
H
T
)
≥
∑
i
k
l
‍
A
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
.
The last term in𝒥
(
H
) is bounded by
(49)
λ
tr
⁡
(
H
T
L
-
H
)
≥
λ
∑
i
k
l
‍
L
k
l
-
H
i
k
′
H
l
i
′
(
1
+
log
⁡
H
i
k
H
l
i
H
i
k
′
H
l
i
′
)
.
Collecting all bounds, we obtainZ
(
H
,
H
′
) as in (41). Obviously, 𝒥
(
H
)
≤
Z
(
H
,
H
′
) and 𝒥
(
H
)
=
Z
(
H
,
H
).
To find the minimum ofZ
(
H
,
H
′
), we take
(50)
∂
Z
(
H
,
H
′
)
∂
H
i
k
=
-
2
B
i
k
+
H
i
k
′
H
i
k
+
2
B
i
k
-
H
i
k
H
i
k
′
+
2
(
H
′
A
+
)
i
k
H
i
k
H
i
k
′
-
2
(
H
′
A
-
)
i
k
H
i
k
′
H
i
k
+
2
λ
(
L
+
H
′
)
i
k
H
i
k
H
i
k
′
-
2
λ
(
L
-
H
′
)
i
k
H
i
k
′
H
i
k
=
0
.
We have(51)
-
2
B
i
k
+
(
H
′
)
i
k
2
+
2
B
i
k
-
H
i
k
2
+
2
(
H
′
A
+
)
i
k
H
i
k
2
-
2
(
H
′
A
-
)
i
k
(
H
′
)
i
k
2
+
2
λ
(
L
+
H
′
)
i
k
H
i
k
2
-
2
λ
(
L
-
H
′
)
i
k
(
H
′
)
i
k
2
=
0
.
Therefore(52)
H
i
k
=
H
i
k
′
B
i
k
+
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
.
The Hessian matrix containing the second derivatives(53)
∂
2
Z
(
H
,
H
′
)
∂
H
i
k
∂
H
j
l
=
δ
i
j
δ
k
l
Y
i
k
is a diagonal matrix with positive entries
(54)
Y
i
k
=
4
[
(
B
+
)
i
k
+
(
H
′
A
-
)
i
k
+
λ
(
L
-
H
′
)
i
k
]
H
i
k
′
H
i
k
2
+
2
B
i
k
-
+
(
H
′
A
+
)
i
k
+
λ
(
L
+
H
′
)
i
k
H
i
k
′
.
Thus,Z
(
H
,
H
′
) is a convex function of H. Therefore, we obtain the global minimum by setting ∂
Z
(
H
,
H
′
)
/
∂
H
i
k
=
0 in (41) and solving for H. Rearranging, we obtain  (21).Theorem 7.
UpdatingV using (21) will monotonically decrease the value of the objective in (13); hence it converges.Proof.
By Lemma2 and Theorem 4, we can get that 𝒥
(
V
t
)
=
Z
(
V
(
t
)
,
V
(
t
)
)
≥
Z
(
V
(
t
+
1
)
,
V
(
t
)
)
≥
𝒥
(
V
t
+
1
), so 𝒥
(
V
) is monotonically decreasing. Since 𝒥
(
V
) is obviously bounded below, we prove this theorem.


## 4. Experimental Results

In this section, we show the performance of the proposed method on face recognition and compare our proposed method with the popular subspace learning algorithms: four unsupervised ones which are principal component analysis [21] (PCA), neighborhood preserving embedding (NPE) [20], local nonnegative matrix factorization (LNMF) [4], and convex nonnegative factorization (CNMF) [13] the one supervised algorithm and which is linear discriminant analysis (LDA) [21]. We use the nearest neighbor (NN) classifier as baseline in original space. We apply different algorithms to obtain new representations for each chosen data set, and then the NN method is applied in the new representation spaces.


### 4.1. Data Preparation

The experiments are used on three data sets. One is Cambridge ORL face database, the other is the Yale database, and the third one is the CMU PIE face database. The important statistics of these data sets are described below.The Yale database contains 165 gray scale images of 15 individuals. All images demonstrate variations in lighting condition (left-light, center-light, right-light), facial expression (normal, happy, sad, sleepy, surprised, and wink), and with/without glasses.The ORL database contains ten different images of each of 40 distinct subjects, thus 400 images in total. For some subjects, the images were taken at different times, varying the lighting, facial expressions (open/closed eyes, smiling/not smiling) and facial details (glasses/no glasses). All the images were taken against a dark homogeneous background with the subjects in an upright, frontal position (with tolerance for some side movement).The CMU PIE face database contains more than 40 000 facial images of 68 people. The images were acquired over different poses, under variable illumination conditions, and with different facial expressions. In our experiment, we choose the images from the frontal pose (C27) and each subject has around 49 images from varying illuminations and facial expressions.In all the experiments, images are preprocessed so that faces are located. Original images are first normalized in scale and orientation such that the two eyes are aligned at the same position. Then the facial areas were cropped into the final images for clustering. Each image is of32
×
32 pixels with 256 gray levels per pixel.


### 4.2. Parameter Settings

For each data set, we randomly divide it into training and testing sets, and evaluate the recognition accuracy on the testing set. In detail, for each individual in the ORL and Yale data sets,we randomly select 2, 3, and 4 images per individual, respectively, for training samples, and the remaining for test samples, while for each individual in the PIE data set, we randomly select 5, 10, and 20 images per individual for training samples. For each partition, we repeated each experiment 20 times and calculated the average recognition accuracy. In general, the recognition rate varies with the dimension of the face subspace. The best result obtained in the optimal subspace and the corresponding dimensionality for each method are shown.For the face recognition experiments, several parameters need to be decided beforehand. For LDA, we use PCA as a first step dimensionality reduction algorithm to avoid the singularity problem. The dimension of the PCA step is fixed asN
-
c and then performs LDA. There are two parameters in our NPCNMF and NPE approach: the number of nearest neighbors k and the regularization parameter λ. Throughout our experiments, we empirically set the number of nearest neighbors k to 5, the value of the regularization parameter λ to 100.Each testing sampley is projected into the linear subspace spanned by the column vectors of the basis matrix U, namely, h
y
≈
W
†
y, where W
† indicates the pseudoinverse of matrix W.


### 4.3. Classification Results

Tables1, 2, and 3 show the evaluation results of all the methods on the three data sets, respectively, where the value in each entry represents the average recognition accuracy of 20 independent trials, and the number in brackets is the corresponding projection dimensionality. These experiments reveal a number of interesting points.(1)
It is clear that the use of dimensionality reduction is beneficial in face recognition. There is a significant increase in performance from using LDA, NPE, NMF, LNMF, and CNMF. However, PCA fails to gain improvement over the baseline. This is because that PCA does not encode the discriminative information.
(2)
The performances of nonnegative algorithms NMF, LNMF, and CNMF are much worse than supervised algorithms LDA, which shows that without considering the labeled data, nonnegative algorithms could not guarantee good discriminating power.
(3)
Our NPCNMF algorithm outperforms all other five methods. The reason lies in the fact that NPCNMF considers the geometrical structure of the data and achieves better performance than the other algorithms. This shows that by leveraging the power of both the parts-based representation and the intrinsic geometrical structure of the data, NPCNMF can learn a better compact representation in the sense of semantic structure.Table 1
Face recognition accuracy on the ORL data set. The number in brackets is the corresponding projection dimensionality.
Method
2 Train
3 Train
4 Train
Baseline
69.32%
77.56%
83.48%
PCA
69.32% (79)
77.56% (118)
83.48% (152)
LDA
72.80% (25)
83.79% (39)
90.13% (39)
NPE
73.19% (36)
84.29% (54)
91.06% (73)
NMF
70.87% (97)
78.98% (81)
84.48% (95)
LNMF
71.73% (178)
81.09% (168)
86.31% (195)
CNMF
72.23% (138)
83.58% (143)
89.56% (111)
NCPNMF
77.31% (143)
86.73% (153)
93.35% (145)Table 2
Face recognition accuracy on the Yale data set. The number in brackets is the corresponding projection dimensionality.
Method
2 Train
3 Train
4 Train
Baseline
46.04%
49.96%
55.62%
PCA
46.04% (29)
49.96% (44)
55.62% (58)
LDA
42.81% (11)
60.33% (14)
68.10% (13)
NPE
48.19% (13)
62.00% (19)
69.00% (73)
NMF
44.11% (112)
49.00% (195)
52.19% (164)
LNMF
44.00% (157)
48.84% (198)
53.57% (197)
CNMF
49.72% (125)
59.50% (168)
65.77% (129)
NPCNMF
63.45% (124)
71.83% (148)
81.38% (153)Table 3
Face recognition accuracy on the PIE data set. The number in brackets is the corresponding projection dimensionality.
Method
5 Train
10 Train
20 Train
Baseline
43.02%
62.90%
83.19%
PCA
42.87% (199)
62.51% (195)
82.84% (200)
LDA
84.39% (67)
90.47% (67)
93.98% (67)
NPE
84.71% (166)
91.48% (200)
94.33% (200)
NMF
78.66% (200)
88.98% (200)
92.52% (200)
LNMF
76.47% (200)
87.91% (200)
92.61% (196)
CNMF
83.72% (176)
90.89% (187)
93.78% (159)
NPCNMF
88.43% (147)
94.86% (158)
98.58% (133)


## 4.1. Data Preparation

The experiments are used on three data sets. One is Cambridge ORL face database, the other is the Yale database, and the third one is the CMU PIE face database. The important statistics of these data sets are described below.The Yale database contains 165 gray scale images of 15 individuals. All images demonstrate variations in lighting condition (left-light, center-light, right-light), facial expression (normal, happy, sad, sleepy, surprised, and wink), and with/without glasses.The ORL database contains ten different images of each of 40 distinct subjects, thus 400 images in total. For some subjects, the images were taken at different times, varying the lighting, facial expressions (open/closed eyes, smiling/not smiling) and facial details (glasses/no glasses). All the images were taken against a dark homogeneous background with the subjects in an upright, frontal position (with tolerance for some side movement).The CMU PIE face database contains more than 40 000 facial images of 68 people. The images were acquired over different poses, under variable illumination conditions, and with different facial expressions. In our experiment, we choose the images from the frontal pose (C27) and each subject has around 49 images from varying illuminations and facial expressions.In all the experiments, images are preprocessed so that faces are located. Original images are first normalized in scale and orientation such that the two eyes are aligned at the same position. Then the facial areas were cropped into the final images for clustering. Each image is of32
×
32 pixels with 256 gray levels per pixel.


## 4.2. Parameter Settings

For each data set, we randomly divide it into training and testing sets, and evaluate the recognition accuracy on the testing set. In detail, for each individual in the ORL and Yale data sets,we randomly select 2, 3, and 4 images per individual, respectively, for training samples, and the remaining for test samples, while for each individual in the PIE data set, we randomly select 5, 10, and 20 images per individual for training samples. For each partition, we repeated each experiment 20 times and calculated the average recognition accuracy. In general, the recognition rate varies with the dimension of the face subspace. The best result obtained in the optimal subspace and the corresponding dimensionality for each method are shown.For the face recognition experiments, several parameters need to be decided beforehand. For LDA, we use PCA as a first step dimensionality reduction algorithm to avoid the singularity problem. The dimension of the PCA step is fixed asN
-
c and then performs LDA. There are two parameters in our NPCNMF and NPE approach: the number of nearest neighbors k and the regularization parameter λ. Throughout our experiments, we empirically set the number of nearest neighbors k to 5, the value of the regularization parameter λ to 100.Each testing sampley is projected into the linear subspace spanned by the column vectors of the basis matrix U, namely, h
y
≈
W
†
y, where W
† indicates the pseudoinverse of matrix W.


## 4.3. Classification Results

Tables1, 2, and 3 show the evaluation results of all the methods on the three data sets, respectively, where the value in each entry represents the average recognition accuracy of 20 independent trials, and the number in brackets is the corresponding projection dimensionality. These experiments reveal a number of interesting points.(1)
It is clear that the use of dimensionality reduction is beneficial in face recognition. There is a significant increase in performance from using LDA, NPE, NMF, LNMF, and CNMF. However, PCA fails to gain improvement over the baseline. This is because that PCA does not encode the discriminative information.
(2)
The performances of nonnegative algorithms NMF, LNMF, and CNMF are much worse than supervised algorithms LDA, which shows that without considering the labeled data, nonnegative algorithms could not guarantee good discriminating power.
(3)
Our NPCNMF algorithm outperforms all other five methods. The reason lies in the fact that NPCNMF considers the geometrical structure of the data and achieves better performance than the other algorithms. This shows that by leveraging the power of both the parts-based representation and the intrinsic geometrical structure of the data, NPCNMF can learn a better compact representation in the sense of semantic structure.Table 1
Face recognition accuracy on the ORL data set. The number in brackets is the corresponding projection dimensionality.
Method
2 Train
3 Train
4 Train
Baseline
69.32%
77.56%
83.48%
PCA
69.32% (79)
77.56% (118)
83.48% (152)
LDA
72.80% (25)
83.79% (39)
90.13% (39)
NPE
73.19% (36)
84.29% (54)
91.06% (73)
NMF
70.87% (97)
78.98% (81)
84.48% (95)
LNMF
71.73% (178)
81.09% (168)
86.31% (195)
CNMF
72.23% (138)
83.58% (143)
89.56% (111)
NCPNMF
77.31% (143)
86.73% (153)
93.35% (145)Table 2
Face recognition accuracy on the Yale data set. The number in brackets is the corresponding projection dimensionality.
Method
2 Train
3 Train
4 Train
Baseline
46.04%
49.96%
55.62%
PCA
46.04% (29)
49.96% (44)
55.62% (58)
LDA
42.81% (11)
60.33% (14)
68.10% (13)
NPE
48.19% (13)
62.00% (19)
69.00% (73)
NMF
44.11% (112)
49.00% (195)
52.19% (164)
LNMF
44.00% (157)
48.84% (198)
53.57% (197)
CNMF
49.72% (125)
59.50% (168)
65.77% (129)
NPCNMF
63.45% (124)
71.83% (148)
81.38% (153)Table 3
Face recognition accuracy on the PIE data set. The number in brackets is the corresponding projection dimensionality.
Method
5 Train
10 Train
20 Train
Baseline
43.02%
62.90%
83.19%
PCA
42.87% (199)
62.51% (195)
82.84% (200)
LDA
84.39% (67)
90.47% (67)
93.98% (67)
NPE
84.71% (166)
91.48% (200)
94.33% (200)
NMF
78.66% (200)
88.98% (200)
92.52% (200)
LNMF
76.47% (200)
87.91% (200)
92.61% (196)
CNMF
83.72% (176)
90.89% (187)
93.78% (159)
NPCNMF
88.43% (147)
94.86% (158)
98.58% (133)


## 5. Conclusion and Future Work

In this paper, we have presented a novel matrix factorization method called NPCNMF for dimensionality reduction, which respects the local geometric structure. As a result, NPCNMF can discriminate power more than the ordinary NMF and CNMF approaches which only consider the Euclidean structure of the data. Experimental results on face datasets show that NPCNMF provides better representation in the sense of semantic structure.Several challenges remain to be investigated in our future work.(1)
A suitable value ofλ is important to our algorithm. It remains unknown how to do model selection theoretically and efficiently.
(2)
NPCNMF is currently limited to the linear projections, and those nonlinear techniques (e.g., kernel tricks) may further boost the algorithmic performance. We will investigate it in our future work.
(3)
Another further research direction is how to extend the current framework for tensor-based nonnegative data decomposition.
(4)
NPCNMF algorithm is iterative and sensitive to the initialization ofW and H. It is unclear how to choose optimal initialization parameters in a principled manner.



---

*Source: 154942-2014-02-25.xml*