# An Inexact Update Method with Double Parameters for Nonnegative Matrix Factorization

Authors: Xiangli Li; Wen Zhang; Xiaoliang Dong; Juanjuan Shi
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Nonnegative matrix factorization (NMF) has been used as a powerful date representation tool in real world, because the nonnegativity of matrices is usually required. In recent years, many new methods are available to solve NMF in addition to multiplicative update algorithm, such as gradient descent algorithms, the active set method, and alternating nonnegative least squares (ANLS). In this paper, we propose an inexact update method, with two parameters, which can ensure that the objective function is always descent before the optimal solution is found. Experiment results show that the proposed method is effective.


---

## Body

## 1. Introduction

Nonnegative matrix factorization (NMF) [1] is not only a well-known matrix decomposition approach but also an utility and efficient feature extraction technique. NMF was first put forward by Lee and Seung. Recently, NMF has been successfully applied in many fields, including face verification [2], text mining [3], gene expression data [4], blind source separation [5], and signal processing [6].A nonnegative matrixA∈R+m×n is decomposed into two low-rank nonnegative matrices W∈R+m×r and H∈R+r×n such that A is approximately equal to WH, denoted by (1)A≈W×H,where r∈N+ and r≪mn/m+n. W and H mean different things to different applications; for example, in blind source separation (BSS), W and H are called mixing matrix and source signal matrix, respectively.In order to decrease the approximation error betweenA and WH, the Euclidean distance-based model is employed in this paper. Namely, NMF can be expressed as the following optimization form: (2)minW,HFW,H=12A-WHF2,s.t.W≥0,H≥0,where ·F is the Frobenius norm and W≥0(H≥0) means that all elements of the matrix W(H) are nonnegative. Clearly, it is difficult to find a global optimal solution (W,H) because the objective function F(W,H) is nonconvex. Therefore, the remaining issue is how to solve the nonconvex problem.In 2001, Lee and Seung tried to find a local optimal solution instead of a global one and proposed multiplicative update algorithm [7]. Multiplicative update algorithm is widely used as an efficient computational method for NMF. The update rules for (2) are given as follows: (3)Hijk+1=HijkWkTVijWkTWkHkij,Wijk+1=WijkVHk+1TijWkHk+1Hk+1Tij,where k represents the iteration count. Later, many new methods were available to solve NMF in addition to multiplicative update algorithm, such as gradient descent algorithms [8, 9], the active set method [10], and alternating nonnegative least squares (ANLS) [11–13].For the sake of solving the minimum problem (2), Hien et al. [14] have proposed a novel algorithm, in which the update rule contains only a parameter and the parameter selection is considered exact in some cases. Compared with multiplicative update algorithm, the novel algorithm [14] has fast convergence speed and small decomposition error.Some inexact technique is widely applied to large scale optimization problems. Based on this idea, we employ an inexact parameter instead of the exact parameter [14]. In the meantime, we add another parameter to accelerate the speed decrease of objective function. Hence, we present an inexact update algorithm with two parameters. The proposed method also updates the elements of matrices W and H one by one. Under some assumptions, the proposed method ensures that the objective function is always descent before the optimal solution is found.Similar to multiplicative update algorithm, the proposed method has many advantages, including its efficient storage, simple calculation, and good results. In multiplicative update algorithm, the descent property is established by means of an auxiliary function. However, the proposed method is descent by the local monotonicity of a quadratic function. And it has a faster convergent speed. Besides, the main idea of ANLS is that two optimization problems can be solved alternately. By contrast, the proposed method is easier to implement.The remainder of this paper is organized as follows: we present the procedure to update an element of the matricesW and H in Section 2. In Section 3, we give an inexact update method with double parameters for NMF and establish the convergent properties. In Section 4, experimental results demonstrate the validity of the method. Finally, we conclude this paper.


## 2. Algorithm for Updating an Element of Matrix

In this section, we will discuss the procedure to update an element of the matricesW and H. In [14], Wij is adjusted by adding a parameter α: (4)W~ij=Wij+α,where(5)α=0,ifq=0,max⁡-qp,-Wij,ifq>0,-qp,ifq<0,p=∑b=1mHjb2,q=∑b=1mWH-Vib×Hjb.Motivated by the above work, we give the following adjustment toWij by two parameters α and β: (6)W~ij=βWij+α,in which β can be viewed as a constant as well as a function. We define the parameter α by a certain value β, where 0<β≤1. Similar to [14], we have that if a≠i, (7)W~Hab=WHab;otherwise,(8)W~Hib=∑p=1rW~ipHpb=∑p≠jWipHpb+W~ijHjb=∑p≠jWipHpb+βWijHjb+αHjb=WHib+β-1WijHjb+αHjb.Next, we deduce (9)fW~,H=12∑a∑bW~Hab-Vab2=12∑a≠i∑bWHab-Vab2+12∑bW~Hib-Vib2=fW,H+gα,where(10)gα=12α+β-1Wij2p+α+β-1Wijq,p=∑bHjb2,q=∑bWH-Vib×Hjb.In order to ensure that the functionf(W~,H) possesses the descent property, we should define α so that g(α) is nonpositive in W~ij=βWij+α≥0. For a given value β, α is expressed as follows: (11)α=1-βWij,ifq=0,max⁡1-βWij-γ,-βWij,ifq>0,1-βWij-γ,ifq<0.Since g(α) is a quadratic function, γ∈(0,-2q/p) or γ∈(-2q/p,0). In other words, the γ value can be “inexact.” If γ=-q/p, specially, g(α) reaches the minimum value. It is obvious that the decline of f(W~,H) is the steepest; at the moment, α is believed to be “exact.”Lemma 1.
IfW~ij does not satisfy the KKT conditions, then f(W~,H)<f(W,H); otherwise, W~=W; this implies f(W~,H)=f(W,H).Proof.
Clearly,Wij≥0, ((WH-V)HT)ij≥0, and Wij×(((WH-V)HT))ij=0 are the KKT conditions of (2). By the definition of q, we have (12)Wij×q=0,Wij≥0,q≥0.
IfW~ij does not satisfy (12), we can obtain g(α)<0 in the condition of q<0 or q>0 and Wij>0. Therefore, f(W~,H)<f(W,H).
Conversely, ifq=0 and α=(1-β)Wij, we have W~=W; if q>0 and Wij=0, it is easy to see that α=0 and W~=W=0. In either case, f(W~,H)=f(W,H).
Similarly, we can deduce the following update rule ofH and Lemma 2:(13)H~ij=β′Hij+α′.Let (14)u=∑aWai2,v=∑aWai×WH-Vaj.For a given value β′, α′ is expressed as follows: (15)α′=1-β′Hij,ifq=0,max⁡1-β′Hij-γ′,-β′Hij,ifq>0,1-β′Hij-γ′,ifq<0,where γ′∈(0,-2v/u) or γ′∈(-2v/u,0).Lemma 2.
IfH~ij does not satisfy the KKT conditions, then f(W,H~)<f(W,H); otherwise, H~=H; this implies f(W,H~)=f(W,H).


## 3. The Proposed Algorithm

Based on the above analysis, in this section, we report our algorithm as follows.Algorithm 3.
( s . 1)
Give the starting pointW1≥0 and H1≥0; set k=1. (s.2)
UpdateWk by using the following formula: Forj=1,…,r and i=1,…,n(16)W~ij=βWij+α. α is computed from (11). (s.3)
UpdateHk by using the following formula: Fori=1,…,r and i=1,…,m(17)H~ij=β′Hij+α′. α ′ is computed from (15). (s.4)
Letk=k+1. Go to (s.2).In order to ensure the monotonous decrease off(W,H), we give the next theorem, which follows directly from Lemmas 1 and 2 and Algorithm 3.Theorem 4.
Suppose that(Wk,Hk) is generated by Algorithm 3; then f(Wk,Hk) is monotonically decreasing; that is, (18)fWk+1,Hk+1<fWk,Hk.Corollary 5.
Sequencef(Wk,Hk) is bounded; then f(Wk,Hk) is a convergence sequence; namely, there exists a positive constant C such that (19)limk→∞⁡fWk,Hk=C,for all k large enough.It is clear to obtain the above corollary from theorem. Next, the another convergence property of Algorithm3 is given.Theorem 6.
Suppose(W-,H-) is a limit point of (Wk,Hk) and (20)∑b=1mH-jb2>0,j=1,…,r,∑a=1nW-ai2>0,i=1,…,r.Then (W-,H-) is the stationary point of problem (2).Since the above theorem corresponds to [14] (Theorem 3.2.2) and the proof is the same as that in [14] (Theorem 3.2.2), we will not prove it here.


## 4. Numerical Experiments

In this section, we give some numerical experiments of Algorithm3 and compare its behavior with the method of [14] (NNMF). All codes of the computer procedures are written in MATLAB and run in a MATLAB 7.10 and are carried out on a PC (CPU 2.13 GHz, 2 G memory) with Windows 7 operation system environment.In experiments, we will compare the following statistics: the number of iterations (Iter), the CPU time in second (time), the speed of convergence to stationary point, and the minimum value of objective function value (fun). Same as [14], the speed of convergence is measured by using the following formula: (21)ΔW,H=δW,HδW+δH,where(22)ΔW,H=∑i=1n∑a=1rmin⁡Wia,WH-VHTia+∑b=1r∑j=1mmin⁡Hbj,WTWH-Vbj.δ(W) is the number of elements of the set (23)min⁡Wia,WH-VHTia≠0,i=1,…,n,a=1,…,r,and δ(H) is the number of elements of the set (24)min⁡Hbj,WTWH-Vbj≠0,j=1,…,m,b=1,…,r.Δ(W,H) is called a normalized KKT residual.In order to avoid the initial point effect on the numerical result, in every experiment, we use 20 initial points, and every initial point(W1,H1) is randomly generated. We list the average values of Iter, Fun, Time, and Δ(W,H), respectively.In Table1, the relevant parameters of Algorithm 3 are specified as follows: (25)β=β′=12,γ=-1.6qp,γ′=-1.6vu.Table 1
A = a b s rand ⁡ n m , r × a b s rand ⁡ n r , n and β=β′=1/2.
(m,n,r)  Iter Fun Time Δ ( W , H ) (50,50,4) Algorithm3 2585 1.239 e - 016 6.989 2.548 e - 009 NNMF 5021 1.523 e - 014 13.104 2.556 e - 009 (100,50,5) Algorithm3 2448 1.808 e - 017 13.037 1.342 e - 009 NNMF 7074 5.450 e - 015 36.535 1.358 e - 009 (100,100,5) Algorithm3 6970 1.689 e - 017 47.876 1.011 e - 009 NNMF 7506 3.889 e - 015 59.702 1.010 e - 009 (200,200,4) Algorithm3 4451 8.486 e - 018 78.766 6.261 e - 010 NNMF 18033 2.282 e - 015 198.386 6.275 e - 010 (200,200,8) Algorithm3 18930 1.860 e - 021 385.494 3.149 e - 010 NNMF 19459 1.015 e - 015 391.797 3.148 e - 010 (200,300,6) Algorithm3 14356 1.150 e - 017 295.138 3.348 e - 010 NNMF 20542 1.335 e - 015 424.479 3.351 e - 010 (500,100,20) Algorithm3 14356 1.150 e - 017 295.138 3.348 e - 010 NNMF 20542 1.335 e - 015 424.479 3.351 e - 010From Table1, we can clearly see the following two points: (1) Algorithm 3 requires fewer iterations and fewer function factorization errors and is less time consuming; (2) for the speed of convergence, Algorithm 3 is not less than NNMF.In Table2, the relevant parameters of Algorithm 3 are specified as follows: (26)β=β′=1,γ=-1.6qp,γ′=-1.6vu.Table 2
A = a b s rand ⁡ n m , r × a b s rand ⁡ n r , n and β=β′=1.
( m , n , r )  Iter Fun Time Δ ( W , H ) ( 50,50,4 ) Algorithm3 2382 2.429 e - 015 6.022 2.526 e - 009 NNMF 4157 1.106 e - 014 13.452 2.546 e - 009 ( 100,50,5 ) Algorithm3 7556 4.610 e - 015 30.280 1.356 e - 009 NNMF 14323 1.624 e - 014 54.663 1.355 e - 009 ( 200,100,10 ) Algorithm3 4572 1.480 e - 016 95.036 3.369 e - 010 NNMF 8171 4.707 e - 015 141.275 3.380 e - 010 ( 200,200,4 ) Algorithm3 9942 7.142 e - 016 107.563 6.277 e - 010 NNMF 15854 2.482 e - 015 162.257 6.279 e - 010 ( 200,300,6 ) Algorithm3 22006 4.955 e - 016 272.971 3.346 e - 010 NNMF 13113 1.396 e - 015 453.043 3.350 e - 010From Table2, we find that Algorithm 3 is still better than NNMF for most of the test problems. Note that Algorithm 3 differs from the method of [14] only in the values of γ and γ′, in which one is “inexact” and another is “exact.” Moreover, from Tables 1 and 2, the value of Fun in Table 1 is more efficient than Table 2 for Algorithm 3.In Table3, the relevant parameters of Algorithm 3 are specified as follows: (27)β=e-k/m,β′=e-k/n,γ=-1.6qp,γ′=-1.6vu.Table 3
A = a b s rand ⁡ n m , r × a b s rand ⁡ n r , n, β=e-k/m, and β′=e-k/n.
(m,n,r)  Iter Fun Time Δ ( W , H ) ( 50,50,4 ) Algorithm3 1635 1.438 e - 015 4.852 2.525 e - 009 NNMF 3896 1.318 e - 014 10.218 2.540 e - 009 ( 100,50,5 ) Algorithm3 2833 6.945 e - 016 13.588 1.350 e - 009 NNMF 6823 7.085 e - 015 31.309 1.354 e - 009 ( 100,100,5 ) Algorithm3 3835 6.004 e - 016 25.241 1.009 e - 009 NNMF 7168 3.711 e - 015 44.616 1.012 e - 009 ( 200,200,4 ) Algorithm3 10703 3.491 e - 016 104.272 6.272 e - 010 NNMF 16664 2.878 e - 015 179.542 6.274 e - 010 ( 200,200,8 ) Algorithm3 15445 1.220 e - 016 207.822 3.139 e - 010 NNMF 20571 1.358 e - 015 399.597 3.153 e - 010 ( 200,300,6 ) Algorithm3 10973 1.930 e - 016 200.539 3.348 e - 010 NNMF 21923 1.907 e - 015 440.734 3.351 e - 010In Table3, β and β′ are regarded as a function, respectively. It is obvious that average numbers of iterations and the minimum value of objective function value for Algorithm 3 are smaller than those for NNMF, CPU time of Algorithm 3 is less than that of NNMF, and the speed of convergence of Algorithm 3 and NNMF is very close.In Figure1, we compare the performance of the proposed method with NNMF, respectively. The convergence curves are given as follows.Figure 1
Performance profiles.From Figure1, it is obvious that the speed of decrease of the objective function value is faster for Algorithm 3. And the numbers of iterations for Algorithm 3 are smaller than that for NNMF.For all test problems, although the two methods have the same speed of convergence to stationary point, the proposed method requires fewer iterations and fewer function factorization errors and is less time consuming. Consequently, the above experiment results show that Algorithm3 performs better than the method of [14] by choosing the parameters β and γ suitably.


## 5. Conclusions

In this paper, we propose an inexact update method, with two parameters, which can ensure that the objective function is always descent before the optimal solution is found. Experiment results show that the proposed method performs better than the method of [14].



---

*Source: 2173914-2016-11-28.xml*