Hypergraph Regularized Lp Smooth Nonnegative Matrix Factorization for Data Representation

1. Introduction

Data representation plays an important role in information retrieval [1], computer vision [2], pattern recognition [3], and other applied fields [4,5]. The dimensions of data matrices are extremely high in these practical applications. The high-dimensional data can cause not only storage difficulties but also possible dimensional curse. Therefore, it is necessary to find an effective and low-dimensional representation of the original high-dimensional data matrix. Matrix factorization is one of the important data representation techniques, and the typical matrix decomposition methods mainly includes: principal component analysis (PCA) [6], linear discriminant analysis (LDA) [8], singular value decomposition (SVD) [7], nonnegative matrix factorization (NMF) [9,10], and so on.

NMF has been gained popularity through the works of Lee and Seung published in Nature [9] and NIPS [10]. It has been widely applied in clustering [11,12,13], face recognition [14,15], text mining [16,17], image processing [18,19,20], hyperspectral unmixing (HU) [21,22], and other fields [23,24,25,26]. Several NMF versions are presented in order to improve data representation capabilities, by introducing with different regularization terms or constraints into the NMF model. For example, by considering the orthogonality of the factor matrices, Ding et al. [12] presented orthogonal nonnegative matrix tri-factorizators (ONMF) approach. By incorporating the graph regularization term into NMF, Cai et al. [2] presented the graph regularized nonnegative matrix factorization (GNMF) , where a simple nearest neighborhood graph is constructed by considering the pairwise geometric relationships between two sample points. However, the models did not take into account the high-order relationships among multiple sample points. Shang et al. [27] presented a graph dual regularization nonnegative matrix factorization (DNMF) approach, which simultaneously considers the intrinsic geometry matrix structures of both the data manifold and the feature manifold. But, DNMF neglects the high-order relationships among multiple sample points or multiple features. To solve the problem of DNMF, Zeng et al. [28] presented the hypergraph regularized nonnegative matrix factorization (HNMF), which incorporates the hypergraph regularization term into NMF and constructs a hypergraph by considering high-order relationships among multiple sample points. But, HNMF is unable to produce a smooth and precise solution because it does not take into account the smoothness of the basis matrix. Recently, Leng et al. [29] proposed the graph regularized $L_p$ smooth nonnegative matrix factorization (GSNMF) by incorporating the graph regularization term and $L_p$ smoothing term into NMF, which considers the intrinsic geometric structures of the sample data and may produce a smooth and more accurate solution of the optimization problem with the addition of graph regularization term and the $L_p$ smoothing constraint. However, in GSNMF, only the pairwise relationships between two sample points are considered, and the high-order relationships among multiple sample points are ignored.

In this paper, by incorporating hypergraph regularization and $L_p$ smoothing constraint terms into the standard NMF, we propose first a hypergraph regularized $L_p$ smooth nonnegative matrix factorization (HGSNMF). The hypergraph regularization term considers the high-order relationships among multiple samples. The $L_p$ smoothing constraint term takes into account the smoothness of the basis matrix, which has been proven it is significant in data representation [30,31,32]. To solve the optimization problem of the HGSNMF model, we offer an effective optimization algorithm using multiplicative update technique and theoretically prove the convergence of the HGSNMF algorithm. Finally, we conduct comprehensive experiments on four data sets to demonstrate the effectiveness of the proposed method.

The main contributions of this work can be summarized as follows:

1) Considering the complex relationships into pairwise ones in simple graph will inevitably lead to loss of the important information. Therefore, we construct the hypergraph regularization term to better discover the hidden semantics and simultaneously capture the underlying intrinsic geometric structure of high-dimensional spacial data samples. In the construction of a hypergraph, each vertex of the hypergraph represents a data sample, and each vertex with its k nearest neighboring samples forms a hyperedge. Each hyperedge represents the similarity relation among a group of samples with higher similarity.

2) We consider the $L_p$ smoothing constraint of the basis matrix, which not only removes the noise of the basis matrix to make it smooth, but also obtains a smooth and more accurate solution to the optimization problem by combining the advantages of isotropic and diffusion anisotropic smoothing.

3) We solve the optimization problem using an efficient iterative technique and conduct comprehensive experiments to empirically analyze our approach on four data sets, the experimental results validate the effectiveness of our proposed method.

The rest of the paper is organized as follows. In Section 2, we introduce some related works including NMF and hypergraph learning. In Section 3, we propose the novel HGSNMF model in details, as well as its updating rules and prove the convergence of the HGSNMF method. In addition, we analysis the complexity of the proposed method. In Section 4, we conduct extensive experiments to validate the proposed method. At last, we conclude this paper in Section 5.

2. Related work

2.1. Nonnegative matrix factorization

Given a nonnegative matrix $\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_n]\in {\mathbb{R}}_{+}^{m\times n}$, each column ${\mathbf{x}}_i\in {\mathbb{R}}_{+}^m(i=1,2,\dots ,n)$ of $\mathbf{X}$ represents a data point. The purpose of NMF is to decompose a nonnegative matrix into two low-rank nonnegative factor matrices $\mathbf{B}\in {\mathbb{R}}_{+}^{m\times r}$ and $\mathbf{C}\in {\mathbb{R}}_{+}^{r\times n}$, the product of which is an approximation of the original matrix. Specifically, the objective function of NMF is

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\min \parallel \mathbf{X}-\mathbf{BC}\parallel }_F^2,\\ \hfill s.t.\mathbf{B}\ge 0,\mathbf{C}\ge 0\end{array}\end{array}$$

(1)

where ${\parallel ·\parallel }_F$ is the Frobenius norm of a matrix, $\mathbf{B}$ is the basis matrix, and $\mathbf{C}$ is the coefficient matrix (also called encoding matrix). Obviously, the objective function is convex for $\mathbf{B}$ or $\mathbf{C}$, but nonconvex for both $\mathbf{B}$ and $\mathbf{C}$. Therefore, it is difficult to find the global minimum. To get a local optimal minimum, Lee et al. [10] proposed iterative multiplicative updating technique to solve the problem (1) as follows:

$$\begin{array}{c}\hfill {\mathbf{B}}_{ik}\leftarrow {\mathbf{B}}_{ik}\frac{{\left({\mathbf{XC}}^{\top }\right)}_{ik}}{{\left({\mathbf{BCC}}^{\top }\right)}_{ik}},\\ \hfill {\mathbf{C}}_{kj}\leftarrow {\mathbf{C}}_{kj}\frac{{\left({\mathbf{B}}^{\top }\mathbf{X}\right)}_{kj}}{{\left({\mathbf{B}}^{\top }\mathbf{BC}\right)}_{kj}},\end{array}$$

where ${\mathbf{B}}^{\top }$ is the transpose of $\mathbf{B}$.

2.2. Hypergraph Learning

Simple graph only considers the pairwise geometric relationship between data samples, the complex internal structure of the data samples could not be efficiently used. To remedy this defect, the hypergraph takes into account the high-order geometric relationships among multiple samples, which can better capture potential geometric information of the data [33]. Thus, hypergraph learning [34,35,36,37,38] is an extension of simple graph learning theory.

The hypergraph $G=(V,E)$ consider the high-order relationships among multiple vertices and is made up of a non-empty vertex set $V=\{v_1,v_2,\dots ,v_n\}$ and a non-empty hypergraph set $E=\{e_1,e_2,\dots ,e_m\}$. Each element $v_i\in V$ is called a vertex, and each element $e_j\in E$ is a subset of V, which is known as a hyperedge of G. G is a hypergraph defined on V if $e_j\ne$, for $j=1,2,\dots ,M$, and $e_1\cup e_2\cup \cdots \cup e_M=V$.

When constructing a hypergraph, we generate the hyperedge by calculating the Euclidean distance between the k nearest neighbors of each vertex. The parameter k is manually set. An incidence matrix $\mathbf{H}\in {\mathbb{R}}_{+}^{\left|V\right|\times \left|E\right|}$ is used to describe the incidence relationship between a vertex and a hyperedge, which is formalized as $\mathbf{H}(v,e)=1$ if $v\in e$, and $\mathbf{H}(v,e)=0$ otherwise.

Figure 1 gives an illustration the simple graph, the hypergraph $\overline{G}=(\overline{V},\overline{E})$ and an incidence matrix. In the undirected simple graph, if a vertex is in the k nearest neighbors of another vertex, the two vertices are joined together by an edge. The hypergraph $\overline{G}=(\overline{V},\overline{E})$ considers the high-order relationships among multiple vertices and is made up of a non-empty vertex set $\overline{V}=\{v_1,v_2,v_3,v_4,v_5,v_6,v_7,v_8\}$ and a non-empty hypergraph set $\overline{E}=\{e_1=\{v_1,v_2,v_4\},e_2=\{v_3,v_4,v_5,v_6\},e_3=\{v_6,v_7,v_8\}\}$. In Figure 1, the solid nodes stand for the vertices, the node sets marked by the solid line segment and the node sets marked by the ellipses denote the hyperedges. Furthermore, each vertex is connected with at least one hypheredge in hypergraph, which is associated with a weight, and each hyperedge can have multiple vertices.

Each hyperedge e can be assigned with a positive number $\omega \left(e\right)$ that presents the weight of the hyperedges. The degree of a vertex v and the degree of a hyperedge can be expressed as $d\left(v\right)=\sum _{e\in E}\omega \left(e\right)\mathbf{H}(v,e)$ and $\delta \left(e\right)=\sum _{v\in V}\mathbf{H}(v,e)$, respectively. According to [34], the unnormalized hypergraph Laplacian matrix can be expressed as

$$\begin{array}{c}\hfill {\mathbf{L}}_{Hyper}={\mathbf{D}}_v-\mathbf{S}\end{array}$$

(2)

where $\mathbf{S}={\mathbf{HWD}}_e^{-1}{\mathbf{H}}^{\top }$, $\mathbf{H}$ is an incidence matrix, $\mathbf{W}$ is a diagonal weight matrix composed of $\omega \left(e\right)$, ${\mathbf{D}}_v$ and ${\mathbf{D}}_e$ denote the diagonal matrices composed of $d\left(v\right)$ and $\delta \left(e\right)$, respectively.

Recently, Wang et al. [21] presented hypergraph-regularized sparse NMF (HGLNMF) for hyperspectral unmixing with incorporates the sparse term into HNMF. HGLNMF considers the sparse of the coefficient matrix, and the hypergraph can simulate the higher-order relationship between multiple pixels by multiple vertices in its hyperedges. Wu et al. [39] presented nonnegative matrix factorization with mixed hypergraph regularization (MHGNMF), by considering higher-order information among the nodes into consideration. Some scholars apply nonnegative matrix factorization to multiple perspectives, Zhang et al. [40] presented semi-supervised multi-view clustering with dual hypergraph regularized partially shared nonnegative matrix factorization(DHPS-NMF). Huang et al. [41] presented diverse deep matrix factorization with hypergraph regularization for multi-view data representation. To some extent, these approaches focus more on the multiple vertices hypergraph, which reflects the higher-order relationships between the multiple vertices, but they ignore the base matrix or its smoothness. To overcome this deficiency, by fully considering the higher-order relationships of multiple vertices and the smoothness of the basis matrix, we propose the following hypergraph regularized $L_p$ smooth nonnegative matrix factorization.

3. Hypergraph Regularized $L_p$ Smooth Nonnegative Matrix Factorization

In this section, we will describe the proposed HGSNMF approach in detail and the iterative updating rules of two factor matrices. Then the convergence of the proposed iterative updating rules is proved. Finally, the calculation cost of this method is given.

We first give the construction of the hypergraph regularization term. Given a nonnegative data matrix $\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_n]\in {\mathbb{R}}_{+}^{m\times n}$, we expect that if two data samples ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ are close, the corresponding encoding vectors ${\mathbf{c}}_i$ and ${\mathbf{c}}_j$ in the low-dimensional space are also close to each other. Constructing the hypergraph for the coefficient matrix, we encode geometrical information by linking each data sample with its k nearest neighbors and denoting their hypergraph connections with heat kernel weight

$$\begin{array}{c}\hfill \omega \left(e_i\right)=\sum _{x_j\in e_i}exp(-\frac{\left|\right|{\mathbf{x}}_i-{\mathbf{x}}_j{\left|\right|}^2}{{\delta }^2}),\end{array}$$

(3)

where $\delta =\frac{1}{kn}\sum _{i=1}^n\sum _{x_j\in e_i}\left|\right|{\mathbf{x}}_i-{\mathbf{x}}_j\left|\right|$ denotes the average distance among all the vertices.

With the above defined weight matrix $\omega \left(e\right)$, the hypergraph regularization of the matrix $\mathbf{C}$ can be caculated by the following optimization problem:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathfrak{R}& =\frac{1}{2}\sum _{e\in E}\sum _{i,j\in e}\frac{\omega \left(e\right)}{\delta \left(e\right)}\left|\right|{\mathbf{c}}_i-{\mathbf{c}}_j{\left|\right|}^2\hfill \\ & =\mathrm{Tr}\left(\mathbf{C}({\mathbf{D}}_v-\mathbf{S}){\mathbf{C}}^{\top }\right)\hfill \\ & =\mathrm{Tr}\left({\mathbf{CL}}_{Hyper}{\mathbf{C}}^{\top }\right),\hfill \end{array}\end{array}$$

where $\mathrm{Tr}(·)$ denotes the trace of a matrix, ${\mathbf{L}}_{Hyper}$ is the hypergraph Laplacian matrix of the hyergraph G and defined by (2).

3.1. The Objective Function

To discover the intrinsic geometric structure information of a data set and produce a smooth and more accurate solution, we propose the HGSNMF method by incorporating hypergraph regularization and the $L_p$ smoothing constraint into NMF. The objective function of our HGSNMF is defined as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}& min\mathcal{O}=\underset{\mathbf{B},\mathbf{C}}{min}\left|\right|\mathbf{X}-{\mathbf{BC}\left|\right|}_F^2+\alpha \mathrm{Tr}\left(\mathbf{C}{\mathbf{L}}_{Hyper}{\mathbf{C}}^{\top }\right)+{2\mu \left|\right|\mathbf{B}\left|\right|}_p^p,\hfill \\ & s.t.\mathbf{B}\ge 0,\mathbf{C}\ge 0,\hfill \end{array}\end{array}$$

(4)

where $\mathbf{B}$ is the basis matrix, ${\left|\right|\mathbf{B}\left|\right|}_p={\left(\sum _{i=1,j=1}^{m,r}{\left({\mathbf{B}}_{ij}\right)}^p\right)}^{\frac{1}{p}},$$0<p\le 2$ and $p\ne 1$, ${\mathbf{B}}_{ij}$ is the ith row and jth column entry of matrix $\mathbf{B}$, and $\mathbf{C}$ is the coefficient matrix, $\alpha$ and $\mu$ are the positive regularization parameters for balancing the reconstruction error in (4). Hypergraph regularization term and $L_p$ smoothing regularization term are presented in the second term and the third term, respectively. The hypergraph regularization term can more effectively discover the hidden semantics and simultaneously capture the underlying intrinsic geometric structure of the high-dimensional space data. The $L_p$ smoothing constraint of the basis matrix not only removes the noise of the basis matrix to make it smooth, but also obtains a smooth and more accurate solution to the optimization problem by combining the advantages of isotropic and diffusion anisotropic smoothing. Then, we will give the detailed derivation of the updating rules, the theoretical proof of convergence and the analysis of the computational complexity of the HGSNMF method, and provide more comparative experiments.

3.2. Optimization method

The objective function $\mathcal{O}$ in (4) is not convex in both $\mathbf{B}$ and $\mathbf{C}$, so it is unrealistic to find the global optimal solution. Thus, we can only obtain the local optimal solution by using iterative method. There are multiplicative update, projective gradient, alternating direction multiplier and dictionary learning algorithm for solving optimization problems. Because the multiplicative update algorithm has the advantages of convergence, and simple operation and small computation, we use multiplicative update algorithm to solve optimization problems. We can turn the objective function into the following unconstrained objective function by using the Lagrange multiplier:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L=& \left|\right|\mathbf{X}-{\mathbf{BC}\left|\right|}_F^2+\alpha \mathrm{Tr}\left(\mathbf{C}{\mathbf{L}}_{Hyper}{\mathbf{C}}^{\top }\right)+{2\mu \left|\right|\mathbf{B}\left|\right|}^p-\mathrm{Tr}\left({\rm Y}{\mathbf{B}}^{\top }\right)-\mathrm{Tr}(\Lambda {\mathbf{C}}^{\top }),\hfill \end{array}\end{array}$$

where ${\rm Y}=\left[{{\rm Y}}_{ik}\right]$, $\Lambda =\left[{\Lambda }_{kj}\right]$, ${{\rm Y}}_{ik}$ and ${\Lambda }_{kj}$ are the Lagrange multiplier for the constrains ${\mathbf{B}}_{ik}\ge 0$ and ${\mathbf{C}}_{kj}\ge 0$, respectively.

Taking the partial derivatives of L with respect to $\mathbf{B}$ and $\mathbf{C}$, respectively, we have that

$$\begin{array}{c}\hfill \frac{\partial L}{\partial \mathbf{B}}=2{\mathbf{BCC}}^{\top }-2{\mathbf{XC}}^{\top }+2\mu p{\mathbf{B}}^{p-1}-{\rm Y},\end{array}$$

$$\begin{array}{c}\hfill \frac{\partial L}{\partial \mathbf{C}}=2{\mathbf{B}}^{\top }\mathbf{BC}-2{\mathbf{B}}^{\top }\mathbf{X}+2\alpha \mathbf{C}{\mathbf{L}}_{Hyper}-\Lambda .\end{array}$$

By using the Karush-Kuhn-Tucker (KKT) conditions $\frac{\partial L}{\partial \mathbf{B}}=0$, $\frac{\partial L}{\partial \mathbf{C}}=0$, ${{\rm Y}}_{ik}·{\mathbf{B}}_{ik}=0$ and ${\Lambda }_{kj}·{\mathbf{C}}_{kj}=0$, we obtain that

$$\begin{array}{c}\hfill {\left({\mathbf{BCC}}^{\top }\right)}_{ik}{\mathbf{B}}_{ik}+\mu p{\mathbf{B}}_{ik}^{p-1}{\mathbf{B}}_{ik}-{\left({\mathbf{XC}}^{\top }\right)}_{ik}{\mathbf{B}}_{ik}=0,\end{array}$$

(5)

$$\begin{array}{c}\hfill \begin{array}{cc}& {\left({\mathbf{B}}^{\top }\mathbf{BC}\right)}_{kj}{\mathbf{C}}_{kj}+\alpha {\left(\mathbf{C}{\mathbf{L}}_{Hyper}\right)}_{kj}{\mathbf{C}}_{kj}-{\left({\mathbf{B}}^{\top }\mathbf{X}\right)}_{kj}{\mathbf{C}}_{kj}=0.\hfill \end{array}\end{array}$$

(6)

According to (5) and (6), we can obtain the following updating rules for $\mathbf{B}$ and $\mathbf{C}$,

$$\begin{array}{c}\hfill {\mathbf{B}}_{ik}\leftarrow {\mathbf{B}}_{ik}\frac{{\left({\mathbf{XC}}^{\top }\right)}_{ik}}{{\left({\mathbf{BCC}}^{\top }\right)}_{ik}+\mu p{\mathbf{B}}_{ik}^{p-1}}\end{array}$$

(7)

and

$$\begin{array}{c}\hfill {\mathbf{C}}_{kj}\leftarrow {\mathbf{C}}_{kj}\frac{{\left({\mathbf{B}}^{\top }\mathbf{X}\right)}_{kj}+\alpha {\left(\mathbf{CS}\right)}_{kj}}{{\left({\mathbf{B}}^{\top }\mathbf{BC}\right)}_{kj}+\alpha {\left(\mathbf{C}{\mathbf{D}}_v\right)}_{kj}},\end{array}$$

(8)

respectively.

3.3. Convergence Analysis

In this subsection, we give the convergence proof of our proposed HGSNMF in (4) using the updating rules (7) and (8). First of all, we introduce some related Definitions and Lemmas.

Definition 1. 

[10]$G(x,x^{\prime })$ is an auxiliary function of $F\left(x\right)$, if $G(x,x^{\prime })$ satisfies the condition

$$G(x,x^{\prime })\ge F\left(x\right),G(x,x)=F\left(x\right).$$

The auxiliary function plays an important role due to the following Lemma.

Lemma 1. 

[10] If G is an auxiliary function of F, then F is nonincreasing under the updating rule

$$\begin{array}{c}\hfill x^{t+1}=arg\underset{x}{min}G(x,x^t).\end{array}$$

(9)

To prove the convergence of HGSNMF under the updating step for $\mathbf{B}$ in (7), we fix the matrix $\mathbf{C}$. For any element ${\mathbf{B}}_{ik}$ in $\mathbf{B}$, we use ${\tilde{F}}_{ik}$ to denote the part of the objective function $\mathcal{O}$, which is only relevant to ${\mathbf{B}}_{ik}$. The first and second derivatives of $\tilde{F}\left({\mathbf{B}}_{ik}\right)$ are given as follows

$$\begin{array}{c}\hfill {\tilde{F}}_{ik}^{\prime }={\left(\frac{\partial \mathcal{O}}{\partial \mathbf{B}}\right)}_{ik}=2{({\mathbf{BCC}}^{\top }-{\mathbf{XC}}^{\top })}_{ik}+2\mu p{\left({\mathbf{B}}_{ik}\right)}^{p-1}\end{array}$$

and

$$\begin{array}{c}\hfill {\tilde{F}}_{ik}^{\prime \prime }={\left(\frac{{\partial }^2\mathcal{O}}{\partial {\mathbf{B}}^2}\right)}_{ik}=2{\left({\mathbf{CC}}^{\top }\right)}_{kk}+2\mu p(p-1){\left({\mathbf{B}}_{ik}\right)}^{p-2},\end{array}$$

respectively.

Lemma 2. 

The function

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \tilde{G}(x,{\mathbf{B}}_{ik}^t)=\tilde{F}\left({\mathbf{B}}_{ik}^t\right)+{\tilde{F}}^{\prime }\left({\mathbf{B}}_{ik}^t\right)(x-{\mathbf{B}}_{ik}^t)+\frac{{\left({\mathbf{BCC}}^{\top }\right)}_{ik}+\mu p{\left({\mathbf{B}}_{ik}^t\right)}^{p-1}}{{\mathbf{B}}_{ik}^t}{(x-{\mathbf{B}}_{ik}^t)}^2\end{array}\end{array}$$

(10)

is an auxiliary function of ${\tilde{F}}_{ik}\left(x\right)$, which is only relevant to ${\mathbf{B}}_{ik}$.

Proof. 

Since from Taylor series expansion, we have that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tilde{F}}_{ik}\left(x\right)=& \tilde{F}\left({\mathbf{B}}_{ik}^t\right)+{\tilde{F}}^{\prime }\left({\mathbf{B}}_{ik}^t\right)(x-{\mathbf{B}}_{ik}^t)+{\left({\mathbf{CC}}^{\top }\right)}_{kk}{(x-{\mathbf{B}}_{ik}^t)}^2\hfill \\ & +\mu p(p-1){\left({\mathbf{B}}_{ik}^t\right)}^{p-2}){(x-{\mathbf{B}}_{ik}^t)}^2.\hfill \end{array}\end{array}$$

(11)

Clearly, $\tilde{G}(x,x)={\tilde{F}}_{ik}\left(x\right)$. By Definition 1, we just need to prove that $\tilde{G}(x,{\mathbf{B}}_{ik}^t)\ge {\tilde{F}}_{ik}\left(x\right).$ By comparing (10) with (11), we can find that $\tilde{G}(x,{\mathbf{B}}_{ik}^t)\ge {\tilde{F}}_{ik}\left(x\right)$ is equivalent to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{{\left(\mathbf{B}{\mathbf{CC}}^{\top }\right)}_{ik}+\mu p{\left({\mathbf{B}}_{ik}^t\right)}^{p-1}}{{\mathbf{B}}_{ik}^t}\ge {\left({\mathbf{CC}}^{\top }\right)}_{kk}+\mu p(p-1){\left({\mathbf{B}}_{ik}^t\right)}^{p-2}.\end{array}\end{array}$$

(12)

Since $\mathbf{B}\ge 0,\mathbf{C}\ge 0,0<p\le 2$, and $p\ne 1$, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\mathbf{B}{\mathbf{CC}}^{\top }\right)}_{ik}=\sum _{l=1,l\ne k}^r{\mathbf{B}}_{il}^t{\left({\mathbf{CC}}^{\top }\right)}_{lk}+{\mathbf{B}}_{ik}^t{\left({\mathbf{CC}}^{\top }\right)}_{kk}\ge {\mathbf{B}}_{ik}^t{\left({\mathbf{CC}}^{\top }\right)}_{kk}\end{array}\end{array}$$

and

$$p{\left({\mathbf{B}}_{ik}^t\right)}^{p-1}\ge p(p-1){\left({\mathbf{B}}_{ik}^t\right)}^{p-1}.$$

So (12) holds, which implies $\tilde{G}(x,{\mathbf{B}}_{ik}^t)\ge {\tilde{F}}_{ik}\left(x\right)$.□

Theorem 1. 

The objective function $\mathcal{O}$ in (4) is nonincreasing under the updating rule of (7).

Proof. 

Replacing $G(x,x^t)$ in (9) by (10), we can obtain the updating rule

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{B}}_{ik}^{t+1}& ={\mathbf{B}}_{ik}^t-{\mathbf{B}}_{ik}^t\frac{{\tilde{F}}^{\prime }\left({\mathbf{B}}_{ik}^t\right)}{2{\left(\mathbf{B}{\mathbf{CC}}^{\top }\right)}_{ik}+2\mu p{\left({\mathbf{B}}_{ik}^t\right)}^{p-1}}\hfill \\ & ={\mathbf{B}}_{ik}^t\frac{{\left({\mathbf{XC}}^{\top }\right)}_{ik}}{{\left(\mathbf{B}{\mathbf{CC}}^{\top }\right)}_{ik}+\mu p{\left({\mathbf{B}}_{ik}^t\right)}^{p-1}}.\hfill \end{array}\end{array}$$

(13)

Since (10) is an auxiliary function of ${\tilde{F}}_{ik}$, ${\tilde{F}}_{ik}$ is nonincreasing under the updating rule of (7). □

Next, we fix the matrix $\mathbf{B}$. For any element ${\mathbf{C}}_{kj}$ in $\mathbf{C}$, we use ${\overline{F}}_{kj}$ to denote the part of the objective function $\mathcal{O}$, which is only relevant to ${\mathbf{C}}_{kj}$. By calculation

$$\begin{array}{c}\hfill {\overline{F}}_{kj}^{\prime }={\left(\frac{\partial \mathcal{O}}{\partial \mathbf{C}}\right)}_{kj}=2{({\mathbf{B}}^{\top }\mathbf{BC}-{\mathbf{B}}^{\top }\mathbf{X}+\alpha {\mathbf{CL}}_{Hyper})}_{kj}\end{array}$$

and

$$\begin{array}{c}\hfill {\overline{F}}_{kj}^{\prime \prime }={\left(\frac{{\partial }^2\mathcal{O}}{\partial {\mathbf{C}}^2}\right)}_{kj}=2{\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}+2\alpha {\left({\mathbf{L}}_{Hyper}\right)}_{jj}.\end{array}$$

Lemma 3. 

The function

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \overline{G}(x,{\mathbf{C}}_{kj}^t)=\overline{F}\left({\mathbf{C}}_{kj}^t\right)+{\overline{F}}^{\prime }\left({\mathbf{C}}_{kj}^t\right)(x-{\mathbf{C}}_{kj}^t)+\frac{{({\mathbf{B}}^T\mathbf{BC}+\alpha \mathbf{C}{\mathbf{D}}_v)}_{kj}}{{\mathbf{C}}_{kj}^t}{(x-{\mathbf{C}}_{kj}^t)}^2\end{array}\end{array}$$

(14)

is an auxiliary function of ${\overline{F}}_{kj}\left(x\right)$, which is only relevant to ${\mathbf{C}}_{kj}$.

Proof. 

Since from Taylor series expansion, we have that

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\overline{F}}_{kj}\left(x\right)=\overline{F}\left({\mathbf{C}}_{kj}^t\right)+{\overline{F}}^{\prime }\left({\mathbf{C}}_{kj}^t\right)(x-{\mathbf{C}}_{kj}^t)+\left({\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}+\alpha {\left({\mathbf{L}}_{Hyper}\right)}_{jj})\right(x-{\mathbf{C}}_{kj}^t{)}^2.\end{array}\end{array}$$

(15)

Clearly, $\overline{G}(x,x)={\overline{F}}_{kj}\left(x\right)$. According to Definition 1, we only need to prove that $\overline{G}(x,{\mathbf{C}}_{kj}^t)\ge {\overline{F}}_{kj}\left(x\right).$ By comparing (14) with (15), we see that $\overline{G}(x,{\mathbf{C}}_{kj}^t)\ge {\overline{F}}_{kj}\left(x\right)$ is equivalent to

$$\begin{array}{c}\hfill \frac{{({\mathbf{B}}^{\top }\mathbf{BC}+\alpha \mathbf{C}{\mathbf{D}}_v)}_{kj}}{{\mathbf{C}}_{kj}^t}\ge {\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}+\alpha {\left({\mathbf{L}}_{Hyper}\right)}_{jj}.\end{array}$$

(16)

Since $\mathbf{B}\ge 0,\mathbf{C}\ge 0$, we have

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left({\mathbf{B}}^{\top }\mathbf{BC}\right)}_{kj}=\sum _{l=1,l\ne j}^n{\mathbf{C}}_{kl}^t{\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}+{\mathbf{C}}_{kj}^t{\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}\ge {\mathbf{C}}_{kj}^t{\left({\mathbf{B}}^{\top }\mathbf{B}\right)}_{kk}\end{array}\end{array}$$

and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left(\mathbf{C}{\mathbf{D}}_v\right)}_{kj}=\sum _{l=1,l\ne j}^n{\mathbf{C}}_{kl}^t{\left({\mathbf{D}}_v\right)}_{lj}+{\mathbf{C}}_{kj}^t{\left({\mathbf{D}}_v\right)}_{jj}\ge {\mathbf{C}}_{kj}^t{({\mathbf{D}}_v-\mathbf{S})}_{jj}={\mathbf{C}}_{kj}^t{\left({\mathbf{L}}_{Hyper}\right)}_{jj}.\end{array}\end{array}$$

So (16) holds, which implies $\overline{G}(x,{\mathbf{C}}_{kj}^t)\ge {\overline{F}}_{kj}\left(x\right)$.□

Theorem 2. 

The objective function $\mathcal{O}$ in (4)is nonincreasing under the updating rule of (8).

Proof. 

Replacing $G(x,x^t)$ in (9) by (14), we can obtain the updating rule

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathbf{C}}_{kj}^{t+1}& ={\mathbf{C}}_{kj}^t-{\mathbf{C}}_{kj}^t\frac{{\overline{F}}^{\prime }\left({\mathbf{C}}_{kj}^t\right)}{2{({\mathbf{B}}^{\top }\mathbf{BC}+\alpha \mathbf{C}{\mathbf{D}}_v)}_{kj}}\hfill \\ & ={\mathbf{C}}_{kj}^t\frac{{({\mathbf{B}}^{\top }\mathbf{X}+\alpha \mathbf{C}\mathbf{B})}_{kj}}{{({\mathbf{B}}^{\top }\mathbf{BC}+\alpha \mathbf{C}{\mathbf{D}}_v)}_{kj}}.\hfill \end{array}\end{array}$$

(17)

Since (14) is an auxiliary function of ${\overline{F}}_{kj}$, ${\overline{F}}_{kj}$ is nonincreasing under the updating rule of (8). □

Similar to NMF, it is known from Theorem 1 and Theorem 2 that the convergence of the model (4) can be guaranteed under the updating rules of (7) and (8).

Now the specific procedure for finding the local optimal $\mathbf{B}$ and $\mathbf{C}$ of HGSNMF is summarized in Algorithm 1.

Algorithm 1 HGSNMF algorithm.

Input: Data matrix $\mathbf{X}=[{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_n]\in {\mathbb{R}}_{+}^{m\times n}$. The number of neighbors k. The algorithm

parameters r, p and regularization parameters $\alpha$, $\mu$. The stopping criterion $ϵ$, and the maximum

number of iterations maxiter. Let ${\Delta }_0=1e-5$.

Output: Factors $\mathbf{B}$ and $\mathbf{C}$ ;

1: Initialize $\mathbf{B}$ and $\mathbf{C}$;

2: Construct the weight matrix $\mathbf{W}$ using (3), and calculate the

matrix ${\mathbf{D}}_v$, $\mathbf{S}$;

3: for$t=1,2,\cdots ,$maxiter do

4:    Update ${\mathbf{B}}^t$ and Update ${\mathbf{C}}^t$ according to (7), (8), respectively.

5:    Compute the objective function value $\mathcal{O}$ of (4) to denote ${\Delta }_t$.

6:    if $\frac{|{\Delta }_t-{\Delta }_{t-1}|}{{\Delta }_{t-1}}<ϵ,$

Break and return $\mathbf{B},\mathbf{C}$.

7:    end if

8: endfor

3.4. Computational Complexity Analysis

In this section, we analyze the computational complexity of the proposed HGSNMF method in comparison to other methods. The fladd, flmlt and fldiv denote a floating-point addition, multiplication and division, respectively. The “O” notation denotes the computational cost. The parameters n, m, r and k denote the number of sample points, features, factors and nearest neighbors to construct an edge or hyperedge, respectively.

According to the updating rules, we calculate the arithmetic operations of each iteration in HGSNMF and summarized the results of the proposed HGSNMF. From Table 1, it shows that the overall costs of our proposed HGSNMF are $O\left(mnr\right)$.

4. Numerical Experimental

In this section, we compare the results of data clustering on four popular data sets to estimate the performance of the proposed HGSNMF method with the related state-of-the-art methods, such as K-means, NMF [9], GNMF [2], HNMF [28], GSNMF [29], HGLNMF [21]. All tests will be done on the computer with a 2-Ghz Intel(R) Core(TM) i5-2500U CPU and 8-GB memory 64-bit MATLAB R2015a in Windows 10. The stopping criterion shows in with $ϵ={10}^{-5}$, and the maximum number of iterations is set to ${10}^4$.

4.1. Data Sets

The clustering performance is evaluated on four widely used data sets, including COIL20 (https://www.cs.columbia.edu/CAVE/software/softlib/coil-20.php), YALE and ORL (https://www.cad.zju.edu.cn/home/dengcai.Data/data.html), and Georgia (https://www.anefian.com/research/face-reco.htm). The important statistical information of the four data sets is listed in Table 2, with more details given as follow. 1) COIL20. The data set contains $32\times 32$ gray scale images of 20 objects viewed at varying angles, and each object has 72 images.

2) ORL. The data set contains $32\times 32$ gray face images of the 40 human subjects taken at different times, different light conditions, different facial expressions, and different facial details, and each object has 10 images.

3) YALE. The data set contains $32\times 32$ gray scale images of the 15 individuals viewed at different facial expression or configuration, and each object has 11 images.

4) Georgia. The data set contains $32\times 32$ gray face images of the 50 people with cluttered background for each object, and each object has 15 colors JPEG images.

4.2. Evaluation Metrics

Two popular evaluation metrics are used to evaluate: the clustering accuracy (ACC) and the normalized mutual information (NMI) [42], which evaluate the clustering performance by comparing the obtained cluster label of each sample with the label provided by the data set. ACC is defined as follows

$$\begin{array}{c}\hfill ACC=\frac{\sum _{i=1}^n\delta (r_i,map\left(q_i\right))}{n},\end{array}$$

where $r_i$ is the correct label provided by the real data set, and $q_i$ is the clustering label obtained by the clustering result, n is the total number of documents, $\delta (x,y)$ is the delta function that equals one if $x=y$ and equals to zero otherwise, and $map(·)$ is a mapping function that maps each cluster label $q_i$ to a given equivalent label from the data set. The Kuhn-Munkers algorithm [43] is used to find the best mapping.

Suppose two clusters C and $C^{\prime }$ are given, then the mutual information metric $MI(C,C^{\prime })$ can be defined as follows

$$\begin{array}{c}\hfill MI(C,C^{\prime })=\sum _{c_i\in C,c_j^{\prime }\in C{}^{\prime }}p(c_i,c_j^{\prime })log\frac{p(c_i,c_j^{\prime })}{p\left(c_i\right)·p\left(c_j^{\prime }\right)},\end{array}$$

where $p\left(c_i\right)$ and $p\left(c_j^{\prime }\right)$ denote the probabilities that a sample arbitrarily selected from the data set belongs to the clusters $c_i$ and $c_j^{\prime }$ , respectively, and $p(c_i,c_j^{\prime })$ denotes the joint possibility that this arbitrarily selected image belongs to the cluster $c_i$ as well as $c_j^{\prime }$ at the same time. The Normalized Mutual Information (NMI) is as defined follows

$$\begin{array}{c}\hfill NMI(C,C^{\prime })=\frac{MI(C,C^{\prime })}{max(H\left(C\right),H\left(C^{\prime }\right))},\end{array}$$

where C is a set of the true labels, and $C^{\prime }$ is a set of clusters obtained from the clustering algorithms. $H\left(C\right)$ and $H\left(C^{\prime }\right)$ are the entropies of C and $C^{\prime }$, respectively.

4.3. Performance Evaluations and Comparisons

To evaluate the performance of our proposed method, we choose K-means, NMF, GNMF, HNMF, and GSNMF as the comparison clustering algorithms.

1) K-means: The K-means performs clustering in the original data, we want to know whether low-dimension data can improve clustering performance on high-dimension data.

2) NMF [9]: The original NMF represented the data by imposing nonnegative constraints on the factor matrices.

3) GNMF [2]: Based on NMF, it constructs local geometric structure of original data space as regularization term.

4) HNMF [28]: It incorporates the hypergraph regularization term into NMF.

5) GSNMF [29]: It incorporates both graph regularization and $L_p$ smoothing constraint into NMF.

6) HGLNMF [21]: It incorporates both hypergraph regularization and $L_{\frac{1}{2}}$ sparse constraint into NMF.

7) HGSNMF: We proposed the hypergraph regularized $L_p$ smooth nonnegative matrix factorization by incorporating hypergraph regularization term and $L_p$ smoothing constraint into NMF.

The general clustering number k is fixed, and we described the experiment as follows.

In NMF, GNMF, HNMF, GSNMF, HGLNMF, and HGSNMF, we initialize two low-rank nonnegative matrix factors using a random strategy in experiments. Next, we set the dimensionality of the low-dimensional space to the number of clusters, and use the classical K-means method to cluster the samples in the new data representation space.

For NMF, GNMF, HNMF, GSNMF, HGLNMF, and HGSNMF, we use the Frobenius norm as the reconstruction error of the objective function. In GNMF and GSNMF, the heat kernel weighting scheme is used to generate the 5 nearest neighbors graph for constructing the weight matrix. In HNMF, HGLNMF, and HGSNMF, the heat kernel weighting scheme is used to generate the 5 nearest neighbors graph for constructing the weight matrix in hypergraph. In GNMF and GSNMF, the graph regularization parameter are all set to 100, in HNMF and HGLNMF, the hypergraph regularization parameter are all set to 100, and in GSNMF, set the parameter $p=1.7$. In HGSNMF, the parameter $\alpha$ and $\beta$ tune from $\{100,500,1000\}$ for COIL20, $\{0.1,10,500,1000\}$ for YALE, $\{1,10,100\}$ for ORL, $\{500,1000\}$ for Georgia; the parameter p tune from $\{0.3,0.5,1.1\}$ for COIL20, $\{0.1,0.5,1.1,1.3\}$ for YALE, $\{0.1,1.2,1.5\}$ for ORL, $\{0.5,0.6\}$ for Georgia, and the best results are reported.

For K-means, we repeat K-means clustering 10 times on original data. For NMF, GNMF, HNMF, GSNMF, HGLNMF, and HGSNMF, we first repeat NMF, GNMF, HNMF, GSNMF, HGLNMF, and HGSNMF 10 times on original data and then repeat K-means clustering 10 times on this low-dimensional reduced data, respectively. We reported the average clustering performance and standard deviation in Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10 and Table 11, where the best results are by bold numbers.

From these experimental results in Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10 and Table 11, we have the following conclusions.

1) The clustering performance of the proposed HGSNMF method on all of the data is better than that of the other algorithms in most cases, which show that the HGSNMF method can find more discriminative information for data. For COIL20, YALE, ORL, and Georgia, the average clustering ACC of the HGSNMF method is more than 1.99%, 1.89%, 0.54%, and 3.37% higher than the second-best method, respectively, and the average clustering NMI of the HGSNMF method is more than 1.84%, 2.26%, 0.6%, and 2.21% higher than the second-best method, respectively,

2)The HGSNMF method is better than the HNMF method in most cases, because the $L_p$ smoothing constraint can combine the merits of isotropic and anisotropic diffusion smooth and more accurate solution.

3)The HGSNMF method is also better than the GSNMF method in most cases, which indicate that the hypergraph regularization term can better discover the underlying geometric information than simple graph regularization term.

For the different cluster numbers on YALE and Geogria data sets, we examine the calculation time based on NMF, GNMF, HNMF, GSNMF, and HGSNMF. In these experiments, we choose the aforementioned identical conditions including parameters and iteration times. From Table 12 and Table 3, NMF has the smallest calculation time since it does not have regularization term. HNMF, GSNMF, and HGSNMF extend GNMF by adding hypergraph regularization, $L_p$ smoothing constraint, and the above two terms respectively, the computation time of HNMF, GSNMF, and HGSNMF is more than GNMF. However, the computation time of HGSNMF is less than HNMF in most cases, which shows the computational advantage of HGSNMF under the $L_p$ smoothing constraint. On the YALE data set, the computation time of HGSNMF is smaller than GSNMF, which indicates that the hypergraph regularization term speeds up the convergence of the proposed HGSNMF method.

4.4. Parameters Selection

There are three parameters in our proposed HGSNMF algorithm: regularization parameters $\alpha ,\mu ,$ and p. When $\alpha ,\mu$ and p are 0, the proposed HGSNMF method reduces to the NMF [9]; when $\mu$ and p are 0, the proposed HGSNMF method reduces to the HNMF [10]. In the experiments, we set $k=5$ for all graph-based and hypergraph-based methods on all data sets. In order to test the effect of each varying parameter, we fix the other varying parameters as subsection 4.3.

Firstly, we adjust the parameter $\alpha$ for GNMF, GSNMF, HNMF, and HGSNMF methods. In HGSNMF, for $k=4$, we set $\mu =100$ and $p=1.1$ in HGSNMF, for $k=8$, we set $\mu =1000$ and $p=0.3$ on COIL20 data set; for $k=11$ and $k=14$, we set $\mu =1000$ and $p=0.5$ on YALE data set; for $k=20$, we set $\mu =1$ and $p=1.5$, for $k=40$, we set $\mu =100$ and $p=1.5$ on ORL data set; for $k=15$, we set $\mu =500$ and $p=0.6$ , for $k=20$ , we set $\mu =1000$ and $p=0.5$ on Geogria data set. Figure 2 demonstrate the accuracy and the normalized mutual information varying with respect to $\alpha$ on four data sets. As can be seen from the Figure 2, HGSNMF performs better than other algorithms in most cases. We can see that the performance of HGSNMF is relatively stable with respect to the parameter $\alpha$ on some data sets.

Next, we adjust the parameter $\mu$, for GSNMF and HGSNMF on four data sets, and set $\alpha =100$ and $p=1.7$ in GSNMF. In HGSNMF, for $k=6$, we set $\alpha =1000$ and $p=1.1$, for $k=12$, we set $\alpha =1000$ and $0.5$ on COIL20 data set; for $k=3$, we set $\alpha =500$ and $p=0.1$, for $k=15$, we set $\alpha =500$, $\alpha =1000$ and $p=0.5$ on YALE data set; for $k=15$, we set $\alpha =10$, $p=0.1$, for $k=25$, we set $\mu =1$ and $p=1.2$ on ORL data set; for $k=5$ and $k=30$, we set $\mu =1000$ and $p=0.6$ on Georgia data set. As can be seen from Figure. Figure 3, HGSNMF performs better than GSNMF in most cases on four data sets, and the performance of HGSNMF is stable with respect to the parameter $\mu$ on some data sets.

Finally, we consider the variation of parameter p. In HGSNMF, for $k=8$, we set $\alpha =\mu =1000$, for $k=10$, we set $\alpha =1000$ and $\mu =500$ on COIL20 data set; for $k=3$, we set $\alpha =0.1$ and $\mu =500$ ; for $k=5$, we set $\alpha =0.1$ and $\mu =1000$ on YALE data set; for $k=5$, we set $\alpha =\mu =10$, for $k=15$, we set $\alpha =10$ and $\mu =1$ on ORL data set; for $k=35$, and $k=45$, we set $\alpha =500$ and $\mu =1000$ on Georgia data set. As show in Figure 4, the performance of HGSNMF is relatively stable and very good with respect to the parameter p varying from 0.1 to 2 on some data sets.

4.5. The Converage Analysis

As described in Section 3, the converge of the proposed HGSNMF method is theoretically proved. In this section, we analyze the convergence of the proposed method through experiments. Figure 5 shows the convergence curves of our HGSNMF method on all four data sets. As can be seen from Figure 5, the objective function is monotonically decreasing and tends to converge after 300 iterations.

5. Conclusion

In this paper, we proposed a hypergraph regularized $L_p$ smooth constrain NMF method for data representation, by incorporating hypergraph regularization term and $L_p$ smoothing constraint into NMF. The hypergraph regularization term can better capture the intrinsic geometry structure of the high dimension space data more comprehensively than simple graph, the $L_p$ smoothing constraint may produce a smooth and more accurate solution to the optimization problem. We present the updating rulers and prove the convergence of our HGSNMF method. Experimental results on four real-world data sets show that our proposed method can achieve better clustering effect than other state-of-the-art methods in most cases.