# Graph Sparse Nonnegative Matrix Factorization Algorithm Based on the Inertial Projection Neural Network

Authors: Xiangguang Dai; Chuandong Li; Biqun Xiang
Journal: Complexity
 (2018)

---

## Abstract

We present a novel method, called graph sparse nonnegative matrix factorization, for dimensionality reduction. The affinity graph and sparse constraint are further taken into consideration in nonnegative matrix factorization and it is shown that the proposed matrix factorization method can respect the intrinsic graph structure and provide the sparse representation. Different from some existing traditional methods, the inertial neural network was developed, which can be used to optimize our proposed matrix factorization problem. By adopting one parameter in the neural network, the global optimal solution can be searched. Finally, simulations on numerical examples and clustering in real-world data illustrate the effectiveness and performance of the proposed method.


---

## Body

## 1. Introduction

Dimensionality reduction plays a fundamental role in image processing, and many researchers have been seeking effective methods to solve this problem. For a given image database, there are many distinct features, whereas the available features are far less enough. Thus, it is of great significance to find useful features with low-dimensionality to represent the original feature space. For this purpose, matrix factorization techniques have attracted great attention in recent decades [1–3], and many different methods have been developed by using different criteria. The most familiar methods include Singular Value Decomposition (SVD) [4], Principal Component Analysis (PCA) [5], and Vector Quantization (VQ) [6]. The main idea of matrix factorization methods is finding several matrices whose product approximates to the original matrix. In dimensionality reduction, the dimension of the decomposed matrices is smaller than that of the original matrix. This gives rise to a low-dimensional compact representation of the original data points, which can facilitate clustering or classification.Among these matrix factorization methods, one of the most used methods is nonnegative matrix factorization (NMF) [3], which requires the decomposed matrices to be nonnegative. The effect of the nonnegative constraint leads NMF to learn a part-based representation of high-dimensional data, and it is applied to so many areas such as signal processing [7], data mining [8, 9], and computer vision [10]. In general, NMF is shown to be effective for unsupervised learning problems, but not applicable to supervised learning problems. To overcome this problem, some researchers [11–13] have presented semi-supervised learning theory to achieve better performance in dimensionality reduction. In the light of locality preserving projection, a graph regularized nonnegative matrix factorization method (GNMF) has been proposed to impose the geometrical information on the data space. The geometrical structure is constructed by a nearest neighbor graph [11]. Based on the idea of label propagation, Liu et al. [13] imposed the label information constraint into nonnegative matrix factorization (CNMF). The idea of CNMF is that the neighboring data points with the same class are supposed to merge together in the low-dimensional representation space.Motivated by previous researches in matrix factorization, in this paper, we propose a novel method, called graph sparse nonnegative matrix factorization (GSNMF), for dimensionality reduction, which can be used for semi-supervised learning problems. In GSNMF, a sparse constraint is imposed on GNMF, and this leads matrix factorization to learn a sparse part-based representation of the original data space. The sparse constraint causes GSNMF to be a nonconvex nonsmooth problem, and traditional optimization algorithms can not be optimized directly. Recently, numerous neural networks have emerged as a powerful tool for optimization problems [14–27]. For some nonconvex problems, an inertial projection neural network (IPNN) [16] has been proposed to search different local optimal solutions by the inertial term. In [17], a shuffled frog leaping algorithm (SFLA) has been developed using the recurrent neural network. Based on the SFLA framework, the global optimal solution can be searched. Moreover, there are many optimization methods for nonconvex nonsmooth problems that use neural networks [22–27].It is worth highlighting some advantages of our proposed method as follows:(i)
Traditional algorithms for GNMF [11] and NMF [3] can easily trap into local optimum solution, and these algorithms are sensitive to initial values, while our proposed algorithm using inertial projection neural network can avoid these problems.(ii)
Our proposed algorithm can be initialized by sparser matrices; however, GNMF and NMF may fail in this case.(iii)
By adopting one parameter in the neural network, GSNMF has the better clustering effect than GNMF and NMF.The rest of the paper is organized as follows. In Section2, some related works to NMF are briefly reviewed; then we introduce GSNMF. Section 3 reviews the inertial projection neural network theory and provides a convergence proof to GSNMF. Section 4 presents numerical examples to demonstrate the validity of our proposed algorithm. Experiments on clustering are given in Section 5. Finally, we present some concluding marks and future work in Section 6.


## 2. Problem Formulation

To find the effective features of high dimensionality data, matrix factorization can be used to learn sets of features to represent data. Given a data matrixY=[Y1,…,Yn]∈Rm×n and an integer k, matrix factorization is to find two matrices A∈Rm×k and S=[S1,…,Sn]∈Rk×n such that(1)Y≈AS.When k≪min⁡m,n, the matrix factorization method can be regarded as a dimensionality reduction method. In image dimensionality reduction, each column of S is a basis vector to capture the original image data and each column of A is the representation with respect to the new basis. The most used method to measure the approximation is the Frobenius norm in the following form:(2)minA,SfA,S=12Y-ASF2=12∑i=1m∑j=1nyij-∑k=1raikskj2.Different matrix factorization methods imposed different constraints on (2) that can solve different practical problems. At present, the most used matrix factorization method is nonnegative matrix factorization (NMF) [3] with nonnegative constraints on A and S. The classic algorithm is summarized as follows:(3)aik=aikYSTikASSTik,skj=skjATYkjATASTkj.Recently, Cai et al. [11] proposed a graph regularized nonnegative matrix factorization method (GNMF) and incorporated the geometrical information into the data space. The goal of GNMF is to find effective basis vectors to represent the intrinsic structure. The research has presented the natural assumption that if two data points yj and yl from Y are close in the intrinsic geometry of the data distribution, new representations of two points sj and sl are also close to each other. For each data point yj, we find its p nearest neighbors and put edges between yj and its neighbors. Edges between each data points yj can be considered as the weight matrix W. If nodes j and l are connected by an edge, then Wjl=1. W can be described by(4)Wjl=1,ifyjandylhave the same label,0,otherwise.The low-dimensional representation ofyj with respect to the new basis is sj. The Euclidean distance is used to measure the dissimilarity between sj and sl by(5)dsj,sl=sj-sl2.With the above analysis, the following term is used to measure the smoothness of the low-dimensional representation:(6)R=12∑j,l=1nsj-sl2Wjl=∑j=1nsjTsjFjj-∑j,l=1nsjTslWjl=tr⁡SFST-tr⁡SWST=tr⁡SLST,where tr⁡(·) denotes the trace of a matrix, Fjj=∑j=1nWjl, and L=F-W. Combining (6) and (2), the new objective function is defined by the Euclidean distance.(7)minA⩾0,S⩾0fA,S=12Y-ASF2+αtr⁡SLST.The algorithm to solve (7) is presented as follows:(8)aik=aikYSTikASSTik,skj=skjATY+SWTkjATAST+SFTkj.Whenα=0 or W=0, GNMF is equivalent to nonnegative matrix factorization. In the representation of the image data, GNMF and NMF only consider the Euclidean structure of image space. However, recent researches have shown that human generated images may from a submanifold of the ambient Euclidean space [28, 29]. In general, the human generated images cannot uniformly fill up the high-dimensional Euclidean space. Therefore, the matrix factorization should respect the intrinsic manifold structure and learn the sparse basis to represent the image data. In the light of sparse coding [30], we impose the sparse constraint on (7) and the optimization problem is transformed into another form:(9)minA⩾0,S⩾0fA,S=12Y-ASF2+αtr⁡SLST+βΦS.Because the optimization problem (9) is nonconvex, a block-coordinate update (BCD) [31] structure is proposed to optimize GSNMF. Given the initial A and S, BCD alternatively solves(10)minS∗fS=12Y-ASF2+αtr⁡SLST+βΦSs.t.S⩾0,(11)minA∗fA=12YT-S∗TATF2s.t.A⩾0until convergence. Since (11) and (10) have a similar form, we only consider how to solve (10); then (11) can be solved accordingly. The problem (10) can be transformed into the following vector form:(12)minHfH=12∑i=1nYi-AHi22+12α∑j,l=1nHj-Hl22Wjl+βH1s.t.H⩾0,where(13)H=H1⋮Hn=S1⋮Sn∈Rnk,Hi∈Rk,Yi∈Rm.It is evident that L1-norm is not differentiable. However, [32] has presented a method to solve it. Supposing H=U-V, problem (12) can be rewritten as follows:(14)minU,VfU,V=12∑i=1nYi-AUi-Vi2+12α∑j,l=1nUj-Vj-Ul-Vl22Wjl+βU-V1s.t.U-V⩾0,where Ui=(Hi)+∈Rk, Vi=(-Hi)+∈Rk, (Hi)+=max⁡Hi,0,i=1,…,n, 1nk=[1,…,1]T, and U and V are, respectively, defined as follows:(15)U=U1⋮Un∈Rnk,V=V1⋮Vn∈Rnk.According to the BCD structure, (14) can be separated into two subproblems. Given the initial U and V, one alternatively solves(16)minUfU∗=12∑i=1nYi+AVi-AUi22+12α∑j,l=1nUj-Ul-Vj-Vl22Wjl+βU1s.t.U⩾V,(17)minVfV∗=12∑i=1nYi-AUi--AVi22+12α∑j,l=1nVj-Vl-Uj-Ul22Wjl+βV1s.t.0≤V≤Uuntil convergence. Since (16) and (17) have a similar form, we only consider how to solve (16); then (17) can be solved accordingly. Equation (16) can be transformed into the following convex quadratic program (CQP):(18)minUfU=12UTCU+DUs.t.U∈Ω,where(19)E=diag⁡ATA,…,ATA∈Rnk×nk,Fij=diag⁡Wij,…,Wij∈Rn×n,F=F11F12⋯F1kF21F22⋯F2k⋮⋮⋮⋮Fk1Fk2⋯Fkk,G=diag⁡F11,F12,…,FkkC=E-αF+αGD=β1nkT-Y1TA+V1TATA,…,YnTA+VnTATA-2αGVT-FVTΩ=U∣U⩾V.According to the above analysis, problem (11) can be also transformed into a convex quadratic program problem. For saving space, we do not provide the derivation process. In the following section, we will introduce IPNN to optimize (18).


## 3. Neural Network Model and Analysis

### 3.1. Inertial Projection Neural Network

To solve problem (18), we establish the following neural network using IPNN [16]:(20)dUdt=X,dXdt=-λX+PΩU-∇fU-U,Ut0=U0∈Ω,where ∇f(U)=CU+DT. Now, we are ready to show the convergence and optimality of (20). For any U∈R2nk and X∈R2nk, we set Q1=(U1T,X1T)T, Q2=(U2T,X2T)T, and R(Qi)=Xi-λXi+PΩUi-∇fUi-Ui; then we will present the following theorems.Theorem 1.
For any initial pointQ0∈Ω with the initial condition Q(0)=Q0, there exists unique continuous solution Q(t)∈Ω, where t∈[t0,+∞).Proof.
Note that∇f(U)=CU+DT is Lipschitz continuous. Let L be the Lipschitz constant. Thus, for any Q1 and Q2∈Ω, we obtain(21)RQ1-RQ2≤X1-X2+λX1-X2+U1-U2+PΩU1-∇fU1-PΩU2-∇fU2≤2+λQ1-Q2+I-CU1-U2≤2+λ+I-CQ1-Q2,where I and L are unit matrix and Lipschitz constant, respectively. Therefore, R(Q) is Lipschitz continuous on Ω. There exists unique solution Q(t) with initial condition Q(t0) by the local existence theorem of solution to ordinary differential equations.Theorem 2.
DefineN(U)=U-PΩ(U-F(U)), if the following two conditions hold.
( 1 ) For any U1 and U2∈Rnk,(22)NU1-NU2,U1-U2⩾cos⁡θ1+I-CNU1-NU2,where θ is the angle between N(U1)-N(U2) and U1-U2.
( 2 ) λ > 1 + I - C / cos ⁡ θ. Then, the solution of model (20) converges to optimal solution set of (18).Proof.
ConsideringU∗∈Ω and the Lyapunov function V(t)=(1/2)Ut-U∗=(1/2)(U(t)-U∗)T(U(t)-U∗), we obtain(23)V¨t+λV˙t=XTX+U-U∗TX˙+λUt-U∗TX=X2+Ut-U∗,X˙+λX=X2+Ut-U∗,PΩU-∇fU-U=X2+Ut-U∗,-NU.Since N(U∗)=0 and Condition (1) holds, (23) can be rewritten as(24)X2=V¨t+λV˙t+Ut-U∗,NU=V¨t+λV˙t+Ut-U∗,NU-NU∗≥V¨t+λV˙t+cos⁡θ1+I-CNU2=V¨t+λV˙t+cos⁡θ1+I-CX˙+λX2=V¨t+λV˙t+cos⁡θ1+I-CX˙2+λcos⁡θ1+I-CdX2dt+λ2cos⁡θ1+I-CX2.Then, (24) can be transformed into(25)V¨t+λV˙t+cos⁡θ1+I-CX˙2+λcos⁡θ1+I-CdX2dt+λ2cos⁡θ1+I-C-1X2≤0,which indicates that the function Zt=V˙t+λVt+λcos⁡θ/1+I-CX2+cos⁡θ/1+I-C∫0tX˙s2ds+λ2cos⁡θ/1+I-C-1∫0tXs2 is monotone nonincreasing. Thus, for any t>0,(26)Zt⩽Z0=V˙0+λV0+cos⁡θ1+I-CX02=U0-U∗TX0+λ2U0-U∗2+λcos⁡θ1+I-CX02.Since λ>1+I-C/cos⁡θ, we obtain λ2cos⁡θ/1+I-C-1>0. Further(27)V˙t+λVt⩽Z0.By multiplying the inequality (27) into eλt, we obtain V˙(t)eλt+λV(t)eλt⩽Z0eλt which implies that(28)dVteλtdt⩽Z0eλt.Integrating (28) from 0 to t, it is obtained that Vt≤V0-1/λZ0e-λt+1/λZ0. Therefore, the trajectory of model (20) is bounded.
SinceV(t) is bounded and inequality (26) holds, we obtain(29)V˙t+cos⁡θ1+I-CXt≤Z0.Thus, (29) can be rewritten as(30)Ut-U∗,Xt+cos⁡θ1+I-CXt2⩽Z0.Since U(t)-U∗ is bounded, (30) indicates that X(t)2 is also bounded. From (25) and (30), one obtains ∫0+∞X(s)2ds<+∞ and ∫0+∞X˙(s)ds<+∞.
Assuming∫0+∞X(s)2ds=Q<+∞, and since ∫0+∞X˙(s)2ds<+∞, it implies that there exists P>0, such that X˙(t)<P for any t∈(0,+∞). Therefore, one obtains(31)∫0+∞Xs2X˙sds≤P∫0+∞Xs2ds<+∞,and it is easy to obtain that(32)∫0+∞Xs2X˙sds≤13limt→+∞⁡Xt3-X03<+∞.
According to the theory of Calculus, the value oflimt→+∞⁡X(t)3 exists. Therefore, the value of limt→+∞⁡X(t) also exists. Since ∫0+∞X(s)2ds<+∞, we obtain limt→+∞⁡X(t)=0. Since V(t) is bounded, we have limt→+∞⁡X(t)=0.
It follows from (20) that(33)limt→+∞⁡X˙+U-PΩU-∇fU=0.If limt→+∞⁡X˙(t)=0, then(34)limt→+∞⁡U-PΩU-∇fU=0,which implies that limt→+∞⁡U(t)=U∗. In the following, we will prove limt→+∞⁡X˙(t)=0.
Definingrε(t)=(1/ε)(X(t+ε)-X(t)), and substituting rε(t) into (20), one obtains(35)r˙εt+λrεt=-1εNUt+ε+1εNUt.Since Condition (1) holds, we have(36)NUt+ε-NUtNUt+ε-NUt≤ε1+I-Ccos⁡θNUt+ε-NUtUt+ε-Utε.Hence, we get(37)NUt+ε-NUt≤ε1+I-Ccos⁡θsups∈t,+∞⁡Xs.
Integrating (35), it gives limt→+∞⁡sup⁡rε(t)=0. Since X˙(t)≤sup⁡rε(t), we obtain limt→+∞⁡X˙(t)=0. Thus, the solution of system (20) converges to the optimal set Ω. The proof is completed.Remark 3.
In Theorem2, Condition (1) should be satisfied. In the following, we discuss the existence of parameter cos⁡θ/1+I-C. For any U1 and U2, we have(38)NU1-NU2=U1-PΩU1-∇fU1-U2+PΩU2-∇fU2≤U1-U2+PΩU1-∇fU1-PΩU2-∇fU2≤U1-U2+I-CU1-U2≤1+I-CU1-U2.Therefore, N(U) is also Lipschitz continuous and(39)NU1-NU22≤1+I-CNU1-NU2U1-U2.Thus, we get(40)NU1-NU22≤1+I-Ccos⁡θNU1-NU2,U1-U2,where θ is the angle between N(U1)-N(U2) and U1-U2. It is easy to obtain(41)NU1-NU2,U1-U2⩾cos⁡θ1+I-CNU1-NU2.


### 3.2. Algorithms

Based on the above analysis, we summarize Algorithms1, 2, and 3 to optimize GSNMF. Firstly, the parameters D, E, F, and G can be derived by Algorithm 1. Secondly, Algorithm 2 applies IPNN to optimize the CQP problem. Thirdly, the optimization problem (9) is divided into two CQP problems which are solved alternatively by Algorithm 3. In the following, we analyse the time complexity of our proposed algorithms. The main cost of our proposed algorithms is spent on the calculation of the gradient ∇f in Algorithm 2. To optimize subproblem (10), the operation in Algorithm 2 is the matrix product EU-UF+UG, which takes O(nk2+n2k). With the cost O(mnk) for calculating DT, the complexity of using Algorithm 2 to optimize the subproblem (10) is (42)Omnk+K1×Onk2+n2k.Similarly, the complexity of using Algorithm 2 to optimize the subproblem (11) is (43)Omnk+K1×Omk2.Hence, the overall cost to solve GSNMF is (44)Omnk+K1×Omk2+nk2+kn2.We summarize the time complexity of one iteration round of GSNMF, GNMF, and NMF in Table 1. At each iteration, there are two O(mnk) operations, same as NMF and GNMF.Table 1
Time complexity of GSNMF, GNMF, and NMF.
Solver Time complexity GSNMF O ( m n k ) + K 1 × O ( m k 2 + n k 2 + k n 2 ) GNMF O ( m n k + m k 2 + n k 2 + k n 2 ) NMF O ( m n k + m k 2 + n k 2 )Algorithm 1:Parameters.
Input: Y, A, V, W, α, β Output: D, E, F, G (1) Calculate the number of rows k and columns n in the matrix V (2) Calculate E with E←ATA (3) Calculate F with F←αW (4) Calculate the diagonal matrix G with Gii←α∑j=1nWij (5) Calculate D with D←β1n×k-(YTA+VTATA)-2α(GVT-FVT)Algorithm 2:IPNN.
Input: U, D, E, F, G, h, λ, Ω, K1 Output: U ∗ Initialization: X, U, i←0 repeat (1) ∇f←EU-UF+UG+DT (2) Ui+1←X∗h+U (3) Xi+1←(-λX+PΩ(U-∇f)-U)∗h+X (4) U=Ui+1 (5) X=Xi+1 (6) i←i+1 until  i≤K1 U∗=UAlgorithm 3:GSNMF.
Input: Y ∈ R m × n, A∈Rm×k, S∈Rk×n, K2 Output: A ∗, S∗ Initialization:  W satifies (6), Ω1=U∣U⩾V, Ω2=V∣0≤V≤U, Ω3=U∣U⩾0, i←0 repeat (1) U=max⁡(S,0), V=max⁡(S,0) (2) [D,E,F,G]←PARAMETER(Y,A,V,W,α,β) (3) U←IPNN(U,D,E,F,G,h,λ,Ω1,K1) (4) [D,E,F,G]←PARAMETER(Y,-A,-U,W,α,β) (5) V←IPNN(V,D,E,F,G,h,λ,Ω2,K1) (6) [D,E,F,G]←PARAMETER(YT,ST,0,0,0,0) (7) A←IPNN(A′,D,E,F,G,h,λ,Ω3,K1) (8) A←AT, i←i+1 until  i≤K2 A∗=A, S∗=S


## 3.1. Inertial Projection Neural Network

To solve problem (18), we establish the following neural network using IPNN [16]:(20)dUdt=X,dXdt=-λX+PΩU-∇fU-U,Ut0=U0∈Ω,where ∇f(U)=CU+DT. Now, we are ready to show the convergence and optimality of (20). For any U∈R2nk and X∈R2nk, we set Q1=(U1T,X1T)T, Q2=(U2T,X2T)T, and R(Qi)=Xi-λXi+PΩUi-∇fUi-Ui; then we will present the following theorems.Theorem 1.
For any initial pointQ0∈Ω with the initial condition Q(0)=Q0, there exists unique continuous solution Q(t)∈Ω, where t∈[t0,+∞).Proof.
Note that∇f(U)=CU+DT is Lipschitz continuous. Let L be the Lipschitz constant. Thus, for any Q1 and Q2∈Ω, we obtain(21)RQ1-RQ2≤X1-X2+λX1-X2+U1-U2+PΩU1-∇fU1-PΩU2-∇fU2≤2+λQ1-Q2+I-CU1-U2≤2+λ+I-CQ1-Q2,where I and L are unit matrix and Lipschitz constant, respectively. Therefore, R(Q) is Lipschitz continuous on Ω. There exists unique solution Q(t) with initial condition Q(t0) by the local existence theorem of solution to ordinary differential equations.Theorem 2.
DefineN(U)=U-PΩ(U-F(U)), if the following two conditions hold.
( 1 ) For any U1 and U2∈Rnk,(22)NU1-NU2,U1-U2⩾cos⁡θ1+I-CNU1-NU2,where θ is the angle between N(U1)-N(U2) and U1-U2.
( 2 ) λ > 1 + I - C / cos ⁡ θ. Then, the solution of model (20) converges to optimal solution set of (18).Proof.
ConsideringU∗∈Ω and the Lyapunov function V(t)=(1/2)Ut-U∗=(1/2)(U(t)-U∗)T(U(t)-U∗), we obtain(23)V¨t+λV˙t=XTX+U-U∗TX˙+λUt-U∗TX=X2+Ut-U∗,X˙+λX=X2+Ut-U∗,PΩU-∇fU-U=X2+Ut-U∗,-NU.Since N(U∗)=0 and Condition (1) holds, (23) can be rewritten as(24)X2=V¨t+λV˙t+Ut-U∗,NU=V¨t+λV˙t+Ut-U∗,NU-NU∗≥V¨t+λV˙t+cos⁡θ1+I-CNU2=V¨t+λV˙t+cos⁡θ1+I-CX˙+λX2=V¨t+λV˙t+cos⁡θ1+I-CX˙2+λcos⁡θ1+I-CdX2dt+λ2cos⁡θ1+I-CX2.Then, (24) can be transformed into(25)V¨t+λV˙t+cos⁡θ1+I-CX˙2+λcos⁡θ1+I-CdX2dt+λ2cos⁡θ1+I-C-1X2≤0,which indicates that the function Zt=V˙t+λVt+λcos⁡θ/1+I-CX2+cos⁡θ/1+I-C∫0tX˙s2ds+λ2cos⁡θ/1+I-C-1∫0tXs2 is monotone nonincreasing. Thus, for any t>0,(26)Zt⩽Z0=V˙0+λV0+cos⁡θ1+I-CX02=U0-U∗TX0+λ2U0-U∗2+λcos⁡θ1+I-CX02.Since λ>1+I-C/cos⁡θ, we obtain λ2cos⁡θ/1+I-C-1>0. Further(27)V˙t+λVt⩽Z0.By multiplying the inequality (27) into eλt, we obtain V˙(t)eλt+λV(t)eλt⩽Z0eλt which implies that(28)dVteλtdt⩽Z0eλt.Integrating (28) from 0 to t, it is obtained that Vt≤V0-1/λZ0e-λt+1/λZ0. Therefore, the trajectory of model (20) is bounded.
SinceV(t) is bounded and inequality (26) holds, we obtain(29)V˙t+cos⁡θ1+I-CXt≤Z0.Thus, (29) can be rewritten as(30)Ut-U∗,Xt+cos⁡θ1+I-CXt2⩽Z0.Since U(t)-U∗ is bounded, (30) indicates that X(t)2 is also bounded. From (25) and (30), one obtains ∫0+∞X(s)2ds<+∞ and ∫0+∞X˙(s)ds<+∞.
Assuming∫0+∞X(s)2ds=Q<+∞, and since ∫0+∞X˙(s)2ds<+∞, it implies that there exists P>0, such that X˙(t)<P for any t∈(0,+∞). Therefore, one obtains(31)∫0+∞Xs2X˙sds≤P∫0+∞Xs2ds<+∞,and it is easy to obtain that(32)∫0+∞Xs2X˙sds≤13limt→+∞⁡Xt3-X03<+∞.
According to the theory of Calculus, the value oflimt→+∞⁡X(t)3 exists. Therefore, the value of limt→+∞⁡X(t) also exists. Since ∫0+∞X(s)2ds<+∞, we obtain limt→+∞⁡X(t)=0. Since V(t) is bounded, we have limt→+∞⁡X(t)=0.
It follows from (20) that(33)limt→+∞⁡X˙+U-PΩU-∇fU=0.If limt→+∞⁡X˙(t)=0, then(34)limt→+∞⁡U-PΩU-∇fU=0,which implies that limt→+∞⁡U(t)=U∗. In the following, we will prove limt→+∞⁡X˙(t)=0.
Definingrε(t)=(1/ε)(X(t+ε)-X(t)), and substituting rε(t) into (20), one obtains(35)r˙εt+λrεt=-1εNUt+ε+1εNUt.Since Condition (1) holds, we have(36)NUt+ε-NUtNUt+ε-NUt≤ε1+I-Ccos⁡θNUt+ε-NUtUt+ε-Utε.Hence, we get(37)NUt+ε-NUt≤ε1+I-Ccos⁡θsups∈t,+∞⁡Xs.
Integrating (35), it gives limt→+∞⁡sup⁡rε(t)=0. Since X˙(t)≤sup⁡rε(t), we obtain limt→+∞⁡X˙(t)=0. Thus, the solution of system (20) converges to the optimal set Ω. The proof is completed.Remark 3.
In Theorem2, Condition (1) should be satisfied. In the following, we discuss the existence of parameter cos⁡θ/1+I-C. For any U1 and U2, we have(38)NU1-NU2=U1-PΩU1-∇fU1-U2+PΩU2-∇fU2≤U1-U2+PΩU1-∇fU1-PΩU2-∇fU2≤U1-U2+I-CU1-U2≤1+I-CU1-U2.Therefore, N(U) is also Lipschitz continuous and(39)NU1-NU22≤1+I-CNU1-NU2U1-U2.Thus, we get(40)NU1-NU22≤1+I-Ccos⁡θNU1-NU2,U1-U2,where θ is the angle between N(U1)-N(U2) and U1-U2. It is easy to obtain(41)NU1-NU2,U1-U2⩾cos⁡θ1+I-CNU1-NU2.


## 3.2. Algorithms

Based on the above analysis, we summarize Algorithms1, 2, and 3 to optimize GSNMF. Firstly, the parameters D, E, F, and G can be derived by Algorithm 1. Secondly, Algorithm 2 applies IPNN to optimize the CQP problem. Thirdly, the optimization problem (9) is divided into two CQP problems which are solved alternatively by Algorithm 3. In the following, we analyse the time complexity of our proposed algorithms. The main cost of our proposed algorithms is spent on the calculation of the gradient ∇f in Algorithm 2. To optimize subproblem (10), the operation in Algorithm 2 is the matrix product EU-UF+UG, which takes O(nk2+n2k). With the cost O(mnk) for calculating DT, the complexity of using Algorithm 2 to optimize the subproblem (10) is (42)Omnk+K1×Onk2+n2k.Similarly, the complexity of using Algorithm 2 to optimize the subproblem (11) is (43)Omnk+K1×Omk2.Hence, the overall cost to solve GSNMF is (44)Omnk+K1×Omk2+nk2+kn2.We summarize the time complexity of one iteration round of GSNMF, GNMF, and NMF in Table 1. At each iteration, there are two O(mnk) operations, same as NMF and GNMF.Table 1
Time complexity of GSNMF, GNMF, and NMF.
Solver Time complexity GSNMF O ( m n k ) + K 1 × O ( m k 2 + n k 2 + k n 2 ) GNMF O ( m n k + m k 2 + n k 2 + k n 2 ) NMF O ( m n k + m k 2 + n k 2 )Algorithm 1:Parameters.
Input: Y, A, V, W, α, β Output: D, E, F, G (1) Calculate the number of rows k and columns n in the matrix V (2) Calculate E with E←ATA (3) Calculate F with F←αW (4) Calculate the diagonal matrix G with Gii←α∑j=1nWij (5) Calculate D with D←β1n×k-(YTA+VTATA)-2α(GVT-FVT)Algorithm 2:IPNN.
Input: U, D, E, F, G, h, λ, Ω, K1 Output: U ∗ Initialization: X, U, i←0 repeat (1) ∇f←EU-UF+UG+DT (2) Ui+1←X∗h+U (3) Xi+1←(-λX+PΩ(U-∇f)-U)∗h+X (4) U=Ui+1 (5) X=Xi+1 (6) i←i+1 until  i≤K1 U∗=UAlgorithm 3:GSNMF.
Input: Y ∈ R m × n, A∈Rm×k, S∈Rk×n, K2 Output: A ∗, S∗ Initialization:  W satifies (6), Ω1=U∣U⩾V, Ω2=V∣0≤V≤U, Ω3=U∣U⩾0, i←0 repeat (1) U=max⁡(S,0), V=max⁡(S,0) (2) [D,E,F,G]←PARAMETER(Y,A,V,W,α,β) (3) U←IPNN(U,D,E,F,G,h,λ,Ω1,K1) (4) [D,E,F,G]←PARAMETER(Y,-A,-U,W,α,β) (5) V←IPNN(V,D,E,F,G,h,λ,Ω2,K1) (6) [D,E,F,G]←PARAMETER(YT,ST,0,0,0,0) (7) A←IPNN(A′,D,E,F,G,h,λ,Ω3,K1) (8) A←AT, i←i+1 until  i≤K2 A∗=A, S∗=S


## 4. Numerical Examples

In this section, we exhibit the global searching ability of GSNMF. By adjusting the inertial termλ in the neural network, different local optimal solutions can be searched. Let m=8, n=8, k=5, 1≤λ≤10, α=0, β=0, K1=5, K2=500, and h=0.2. In order to ensure the validity of this experiment, we provide the initial Y, A, and S in Tables 2, 3, and 4, respectively. Table 5 shows the comparison between GSNMF and NMF. To investigate whether GSNMF can converge, the convergence curve is depicted in Figure 1 with λ=1.Table 2
Y.
4.345 1.79 2.92 0.852 6.036 6.929 2.803 2.445 0.513 1.164 8.54 8.645 0.879 12.309 4.731 1.326 4.277 3.665 5.619 2.1 6.048 6.164 1.077 4.167 5.658 5.788 0.887 4.599 7.653 2.342 7.534 7.104 2.282 0.3 12.657 3.248 4.33 1.508 1.819 1.26 1.856 7.669 2.113 4.084 7.616 0.287 2.865 1.329 0.589 6.886 2.388 0.446 7.759 2.456 0.486 7.189 0.214 16.156 3.684 1.973 0.25 3.841 2.825 1.482Table 3
A.
0.797 0.394 0.001 1.029 0.182 0.3 0.633 0.28 1.14 0.727 0.81 0.078 1.589 1.764 1.152 0.062 0.408 0.265 0.78 0.716 0.073 0.389 0.005 1.768 0.166 1.921 1.45 1.474 0.673 0.016 0.346 1.896 0.284 1.135 0.087 1.849 0.058 0.445 1.574 1.745Table 4
S.
0.071 0.242 0.413 1.101 1.345 0.967 0.923 0.533 1.319 1.305 0.329 0.99 0.309 1.173 0.734 1.309 1.377 0.796 0.913 0.377 0.504 1.322 0.378 0.682 0.403 0.608 0.224 0.163 0.612 0.162 0.724 0.849 0.929 1.402 0.773 1.727 0.039 0.187 0.446 0.132Table 5
Comparisons between GSNMF using differentλ and NMF.
GNMF NMF Objective values λ = 1, 20.5172 20.5172 λ = 3, 20.7534 λ = 6, 21.6804Figure 1
Iteration number versus objective value for GSNMF whenλ=1.


## 5. Application in Image Clustering

### 5.1. Databases

To examine the clustering performance of GSNMF, we present the experiment in two databases including IRIS and COIL20. Their details are presented in the following (see also Table6).Table 6
Two data sets.
Dataset size(n) Dimensionality(m) Classes(k) IRIS 150 4 3 COIL20 1440 256 20(  1) IRIS. It includes 150 instances with 4 features. There are 3 classes including Versicolour, Setosa, and Virginca. Each class includes 50 instances.(  2) COIL20. This data set is an image library which contains 1440 instances with 16 × 16 gray scale features. There are 20 different classes and each class contains 72 instances.


### 5.2. Compared Methods

We present the clustering performance on two databases using GSNMF and GNMF. There are two metrics including accuracy and normalized mutual information [33] to evaluate the clustering performance. To reveal the effect of the sparse constraint, different cardinalities (the number of zero entries) of the initial S are considered to evaluate the clustering performance. Because NMF is a nonconvex problem, different initial A and S may lead to different local optimal solutions. For the fair comparison, we try 10 different initial A and S and report the average results. We compare different methods in two cases. Firstly, A and S are randomly generated between 0 and 1. Secondly, A and S are randomly generated between −1 and 1.


### 5.3. Compared Results

Tables7 and 8 present the clustering results on IRIS in two cases. The parameters for IRIS are h=0.02, K1=3, K2=200, α=1, β=1e-4, and 0≤λ≤10. The clustering performance on COIL20 is shown in Tables 9 and 10 in two cases. The parameters for COIL20 are h=0.02, K1=20, K2=200, α=1, β=1e-4, and 0≤λ≤10. These tables reveal some interesting points:(i)
When the sparse constraint is imposed on GNMF, the clustering performance of GSNMF is better than NMF and GNMF.(ii)
When the initialA and S have some negative entries, NMF and GNMF fail to cluster. However, GSNMF is not affected in this case.Table 7
Clustering performance on IRIS in case1.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 98.27 91.60 83.87 93.09 80.66 70.58 10 98.47 77.93 69.47 93.97 55.70 48.48 20 98.27 71.93 68.67 93.09 42.23 38.86 30 98.33 61.93 62.87 93.39 28.90 28.82 40 98.27 56.93 55.80 93.15 18.89 17.59 50 98.33 53.40 53.27 93.39 13.35 12.96Table 8
Clustering performance on IRIS in case2.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 98.60 34.00 34.00 94.56 1.035 1.035 10 98.13 34.00 34.00 93.39 1.035 1.035 20 98.20 34.00 34.00 93.02 1.035 1.035 30 98.00 34.00 34.00 92.27 1.035 1.035 40 98.13 34.00 34.00 92.92 1.035 1.035 50 98.00 34.00 34.00 92.33 1.035 1.035Table 9
Clustering performance on COIL20 in case1.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 68.92 72.92 64.54 77.52 84.88 74.18 10 68.24 63.37 57.63 77.05 71.78 65.49 20 68.15 53.40 48.39 77.31 61.90 56.26 30 68.38 45.81 43.91 77.39 54.34 51.60 40 68.06 41.03 38.28 77.31 49.45 46.26 50 68.67 35.25 32.64 77.61 42.92 41.42Table 10
Clustering performance on COIL20 in case2.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 68.22 5.07 5.07 77.05 1.38 1.035 10 68.16 5.07 5.07 77.36 1.38 1.035 20 68.78 5.07 5.07 77.11 1.38 1.035 30 68.12 5.07 5.07 77.08 1.38 1.035 40 68.38 5.07 5.07 77.33 1.38 1.035 50 68.25 5.07 5.07 77.67 1.38 1.035


### 5.4. Parameters Selection

Our GSNMF has one essential parameter: the inertia termλ. Figures 2 and 3 depict the average performance of GSNMF with different λ.Figure 2
The performance of GSNMF versus differentλ on IRIS data set.
(a)
Accuracy (b)
Normalized mutual informationFigure 3
The performance of GSNMF versus differentλ on COIL20 data set.
(a)
Accuracy (b)
Normalized mutual information


### 5.5. Convergence Study

According to the BCD and IPNN theory, the method for optimizing GSNMF is proved to be convergent. Here we investigate whether this method can converge to a stationary point. Figure4 depicts the convergence curves of GSNMF on two data sets. For each figure, the x-axis is the iteration number and the y-axis denotes the objective value.Figure 4
Iteration number versus objective value (100 : 1 scale) for GSNMF whenλ=5 on IRIS and COIL20 data sets.
(a)
IRIS (b)
COIL20


## 5.1. Databases

To examine the clustering performance of GSNMF, we present the experiment in two databases including IRIS and COIL20. Their details are presented in the following (see also Table6).Table 6
Two data sets.
Dataset size(n) Dimensionality(m) Classes(k) IRIS 150 4 3 COIL20 1440 256 20(  1) IRIS. It includes 150 instances with 4 features. There are 3 classes including Versicolour, Setosa, and Virginca. Each class includes 50 instances.(  2) COIL20. This data set is an image library which contains 1440 instances with 16 × 16 gray scale features. There are 20 different classes and each class contains 72 instances.


## 5.2. Compared Methods

We present the clustering performance on two databases using GSNMF and GNMF. There are two metrics including accuracy and normalized mutual information [33] to evaluate the clustering performance. To reveal the effect of the sparse constraint, different cardinalities (the number of zero entries) of the initial S are considered to evaluate the clustering performance. Because NMF is a nonconvex problem, different initial A and S may lead to different local optimal solutions. For the fair comparison, we try 10 different initial A and S and report the average results. We compare different methods in two cases. Firstly, A and S are randomly generated between 0 and 1. Secondly, A and S are randomly generated between −1 and 1.


## 5.3. Compared Results

Tables7 and 8 present the clustering results on IRIS in two cases. The parameters for IRIS are h=0.02, K1=3, K2=200, α=1, β=1e-4, and 0≤λ≤10. The clustering performance on COIL20 is shown in Tables 9 and 10 in two cases. The parameters for COIL20 are h=0.02, K1=20, K2=200, α=1, β=1e-4, and 0≤λ≤10. These tables reveal some interesting points:(i)
When the sparse constraint is imposed on GNMF, the clustering performance of GSNMF is better than NMF and GNMF.(ii)
When the initialA and S have some negative entries, NMF and GNMF fail to cluster. However, GSNMF is not affected in this case.Table 7
Clustering performance on IRIS in case1.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 98.27 91.60 83.87 93.09 80.66 70.58 10 98.47 77.93 69.47 93.97 55.70 48.48 20 98.27 71.93 68.67 93.09 42.23 38.86 30 98.33 61.93 62.87 93.39 28.90 28.82 40 98.27 56.93 55.80 93.15 18.89 17.59 50 98.33 53.40 53.27 93.39 13.35 12.96Table 8
Clustering performance on IRIS in case2.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 98.60 34.00 34.00 94.56 1.035 1.035 10 98.13 34.00 34.00 93.39 1.035 1.035 20 98.20 34.00 34.00 93.02 1.035 1.035 30 98.00 34.00 34.00 92.27 1.035 1.035 40 98.13 34.00 34.00 92.92 1.035 1.035 50 98.00 34.00 34.00 92.33 1.035 1.035Table 9
Clustering performance on COIL20 in case1.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 68.92 72.92 64.54 77.52 84.88 74.18 10 68.24 63.37 57.63 77.05 71.78 65.49 20 68.15 53.40 48.39 77.31 61.90 56.26 30 68.38 45.81 43.91 77.39 54.34 51.60 40 68.06 41.03 38.28 77.31 49.45 46.26 50 68.67 35.25 32.64 77.61 42.92 41.42Table 10
Clustering performance on COIL20 in case2.
Cardinality (%) Accuracy (%) Normalized mutual information (%) SGNMF GNMF NMF SGNMF GNMF NMF 0 68.22 5.07 5.07 77.05 1.38 1.035 10 68.16 5.07 5.07 77.36 1.38 1.035 20 68.78 5.07 5.07 77.11 1.38 1.035 30 68.12 5.07 5.07 77.08 1.38 1.035 40 68.38 5.07 5.07 77.33 1.38 1.035 50 68.25 5.07 5.07 77.67 1.38 1.035


## 5.4. Parameters Selection

Our GSNMF has one essential parameter: the inertia termλ. Figures 2 and 3 depict the average performance of GSNMF with different λ.Figure 2
The performance of GSNMF versus differentλ on IRIS data set.
(a)
Accuracy (b)
Normalized mutual informationFigure 3
The performance of GSNMF versus differentλ on COIL20 data set.
(a)
Accuracy (b)
Normalized mutual information


## 5.5. Convergence Study

According to the BCD and IPNN theory, the method for optimizing GSNMF is proved to be convergent. Here we investigate whether this method can converge to a stationary point. Figure4 depicts the convergence curves of GSNMF on two data sets. For each figure, the x-axis is the iteration number and the y-axis denotes the objective value.Figure 4
Iteration number versus objective value (100 : 1 scale) for GSNMF whenλ=5 on IRIS and COIL20 data sets.
(a)
IRIS (b)
COIL20


## 6. Conclusion and Future Work

We propose a dimensionality reduction method, which can be solved by the inertial projection neural network. According to the experiments, three advantages are presented. Firstly, different local solutions can be achieved with different inertial termsλ. Secondly, the clustering performance cannot be affected by the negative initial values. However, GNMF and NMF have poor performance in clustering with negative initial values. Thirdly, if the initial values are sparse, our proposed method performs better than GNMF and NMF in the clustering.Several topics remain to be discussed in our future work:(i)
There is a parameterλ which searches the global optimal solution of GSNMF. Thus, a suitable value of λ is critical to our algorithm. It remains unclear how to select λ theoretically.(ii)
h is a step length to decide the convergence rate in Algorithm 2. If it is assigned a small value, slow convergence makes a bad clustering performance. Thus, an adaptive step length should be considered.



---

*Source: 2743678-2018-03-19.xml*