# A Symmetric Rank-One Quasi-Newton Method for Nonnegative Matrix Factorization

Authors: Shu-Zhen Lai; Hou-Biao Li; Zu-Tao Zhang
Journal: ISRN Applied Mathematics
 (2014)

---

## Abstract

As is well known, the nonnegative matrix factorization (NMF) is a dimension reduction method that has been widely used in image processing, text compressing, signal processing, and so forth. In this paper, an algorithm on nonnegative matrix approximation is proposed. This method is mainly based on  a relaxed active set and the quasi-Newton type algorithm, by using the symmetric rank-one and negative curvature direction technologies to approximate the Hessian matrix. The method improves some recent results. In addition, some numerical experiments are presented in the synthetic data, imaging processing, and text clustering. By comparing with the other six nonnegative matrix approximation methods, this method is more robust in almost all cases.


---

## Body

## 1. Introduction

An NMF problem is to decompose a nonnegative matrixV∈ℝn×m into two nonnegative matrix W∈ℝn×k and H∈ℝk×m, such that WH approximate to V as well as possible. To measure the distance between WH and V, there are many methods such as the Kullback Leibler divergence, and Bregman divergence, and Frobenius divergence. Due to the favorable property of the Frobenius divergence, many methods are presented based on it. The Frobenius divergence may be described as follows:
(1)min⁡(W∈ℝ(n×k),H∈ℝ(k×m))f(W,H)=12∥V-WH∥F2,
subject to Wij≥0, Hij≥0 for all i and j.In the last decade, numerous methods have been proposed to deal with the NMF problem in (1). Most of these methods can be classified into two classes, that is, alternating one-step gradient descent and alternating least squares.(i)
The alternating one-step gradient is to alternateW and H with one step. The most well known is Lee's multiplicative update algorithm [1], which alternates W and H by the following rules: suppose we have obtained the lth matrix Wl and Hl, then
(2)Wial+1=Wial(V(Hl)T)ia(WlHl(Hl)T)ia,∀i,a,Hbjl+1=Hbjl((Wl+1)TV)bj((Wl+1)TWl+1Hk)bj,∀b,j.(ii)
For the frame work of the alternating least squares, the procedure is as follows [2]:
(1)
initializeH∈ℝk×m with a nonnegative matrix;(2)
solve the following problems repeatedly until a stopping criterion is satisfied:(3)min⁡W≥0g(W)=12∥HTWT-VT∥F2,
where H is fixed, and
(4)min⁡H≥0g(H)=12∥V-WH∥F2,
where W is fixed.Obviously, the above two methods both satisfy that(5)f(Wl+1,Hl)≤f(Wl,Hl),f(Wl+1,Hl+1)≤f(Wl+1,Hl).
Equation (5) insures that the above numerical methods converge to the solution of (1).Since (3) and (4) may be regarded as the subproblems of the NMF problem (1). Next we mainly focus on how to solve (4), while (3) is similar to (4), due to the symmetric property of Euclidear distance. There exist some methods to solve (4), such as the projected Newton method [3], quasi-Newton method [4–6], the active set method [7–9], and so forth.Note that the Frobenius norm of a matrix is just the sum of the Eucldean distance over columns (or rows). Then solving (4) in fact is equal to solving a series of the following nonnegative least squares (NNLS) problems:
(6)min⁡h1≥0g(h1)=12∥Wh1-v1∥22,min⁡h2≥0g(h2)=12∥Wh2-v2∥22,…,min⁡hk≥0g(hk)=12∥Whk-vk∥22,
where hi,vi, i=1,2,…,m are the columns of H and V, respectively. But, for convenience, we write the least squares by omitting their subscripts. For example,
(7)min⁡g(h)=12∥Wh-v∥22,
subject to h≥0 for all i.As is well known, ifx is the optimal solution of (7), then x satisfies the following equation:
(8)r=WTWh-WTv,r≥0,x≥0,xTv=0,
which is called the KKT optimal condition [10]. This tells us that one can try to find a descend method whose solutions will satisfy the condition (8) after finite iterations.The remained part of this paper is organized as follows. In the second section we will propose our algorithm for a single right hand vector nonnegative least squares. In addition, the symmetric rank-one quasi-Newton method for NMF will be discussed in the third section. In the fourth part we present the numerical experiments on both synthetic data and real world data. The last part is to be the conclusion part, in this section we will forecast the future work on NMF.


## 2. A New Method for Nonnegative Least Squares Problems

The nonnegative least squares can be regarded as a quadratic programming with box constrains, whose upper constrains can be considered as infinite. Recently, Lin proposed the projected algorithm for nonnegative linear square (NNLS) in [3]. Later in [5], Kim et al. proposed the fast projection-based methods (fnmae) for nonnegative least square programming based on the quasi-Newton and active set method. A few years later, Gong and Zhang [11] improved the fnmae algorithm by using the symmetric property of the Hessian matrix of the object function f(W,H). In this paper, we will use the symmetric rank-one quasi-Newton line search technology approximate the Hessian matrix. In addition, we will combine the block active method with the quasi-Newton method and then use the symmetric rank-one technology to modify the BFGS in Kim et al. [5].


### 2.1. The Frame of Solving NNLS

For the NNLS problem (7), our method is iterative and in each iteration we partition the variables into two parts, namely, the inactive set and the active set. Suppose the current iteration is l, then we are going to compute the (l+1)th iteration.Based on the idea of fixed set in [5] and the theory in [12], we denote an inactive set as follows:
(9)I+l={i∣0≤hil≤ϵl,[∇g(h)]i>0},
where ϵl=min⁡{ϵ,∥hl-∇g(hl)∥F2}, ϵ is a small positive scalar, and [∇g(h)]i denotes the ith element of the gradient ∇g(h) of the object function g(h).For convenience, in some places, we will slightly abuse notations and say thathil∈I+l whenever i∈I+l.Denote the active variables and the inactive variables at the current iteration byξl and ηl, respectively. Assuming that the hl and ∇g(hl) are partitioned as follows:
(10)hl=[ξlηl],
where ξl∉I+l and ηl∈I+l. That is to say ξl is the active part of the current iteration. Then we compute the projection of ξ by the following equation:
(11)ξ=P[ξl-αD¯l∇g(ξl)],
where α≥0 is the iteration step-size, D¯l is a gradient scaling matrix, and P[a] is the positive projection, that is,
(12)P[a]={a,a>00,a≤0.
Next we update the current result hl by the following rules:
(13)hl+1=[ξηl]=[P[ξl-αD¯l∇g(ξl)]0],
where the right equation uses the fact that ηl is fixed to zeros, it can be comprehended as that this part satisfies the KKT optimal condition (8), then in the next iteration we need not to update this part. After updating the hl+1, we can compute the ∇g(hl+1) and update the inactive set I+l+1 to obtain the next iterative result. In the whole iteration we will adhere to the iteration rule until the result satisfies the stop criterions.


### 2.2. The Approximate toD¯l

As the size ofξl and ηl changes at each iteration, the computation of the matrix D¯l is not an easy task. Note that the curvature information of ξl is received from the curvature information of hl. Due to this fact, the gradient matrix D¯l can be obtained by taking the proper submatrix of the matrix D¯ to make this task easier, where D¯ is the matrix that carries the curvature information of the vector hl. In realization of the method we can try to eliminate the curvature information of ηl. Then the update rule of (13) can be regarded as follows:
(14)hl+1=P[hl-αD¯∇g(h)].
In [5], Kim et al. used the BFGS update method to approximate the Hessian matrix of the object function g(h). The BFGS method is well established and only uses the gradient information of the object function. But BFGS needs the time consumption. Many researchers have experimentally observed that the symmetric rank-one (SR1) rule performs better than BFGS quasi-Newton update rules [13]. In this paper, we will use the SR1 rule to approximate the Hessian matrix of g(h) to improve the the converge speed.Suppose thatHl is the current approximation of the Hessian matrix, then we apply the following SR1 rule to approximate the next Hessian matrix:
(15)Hl+1=Hl+(wl-Hlvl)(w-Hlvl)T(w-Hlvl)Twl,
where wl=∇g(hl+1)-∇g(hl) and vl=hl+1-hl. Let Dl be the inverse of the Hl. By using the Sherman-Morrison-Woodbury formula [14], we can obtain the inverse matrix Dl+1 of Hl+1 as follows:(16)Dl+1=Dl+(vl-Dlwl)(vl-Dlwl)T(vl-Dlwl)Twl.
Thus, we complete the approximate to D¯l in the iterative processes.


### 2.3. Line Search Strategy for Step-Size

For (13), α>0, we next consider how to find an α with the largest function reduction. For this purpose, we may choose the searching rule as follows:(17)α=argmin⁡α≥0g(P[hl-α∇g(hl)]).
To sum up Sections 2.1–2.3, we obtain the following symmetric rank-one NNLS algorithm (see Algorithm 1).Algorithm 1:The SR1 algorithm for NNLS.
Step  1. start data
Chooseh0≥0, and the maximum iterative number maxiter,
Setl=0;
Step  2.
Compute the inactive set by using the rule (9);
Step  3. (update the current solution)
Renew current solution by (13), where using (16) to approximate
the inverse of Hessian matrix and (17) to compute the stepsize;
Step  4.
ifthe current result reaches the stop criterion then stop, the current
solution is the final solution;
else go to Step  2 and Step  3.


## 2.1. The Frame of Solving NNLS

For the NNLS problem (7), our method is iterative and in each iteration we partition the variables into two parts, namely, the inactive set and the active set. Suppose the current iteration is l, then we are going to compute the (l+1)th iteration.Based on the idea of fixed set in [5] and the theory in [12], we denote an inactive set as follows:
(9)I+l={i∣0≤hil≤ϵl,[∇g(h)]i>0},
where ϵl=min⁡{ϵ,∥hl-∇g(hl)∥F2}, ϵ is a small positive scalar, and [∇g(h)]i denotes the ith element of the gradient ∇g(h) of the object function g(h).For convenience, in some places, we will slightly abuse notations and say thathil∈I+l whenever i∈I+l.Denote the active variables and the inactive variables at the current iteration byξl and ηl, respectively. Assuming that the hl and ∇g(hl) are partitioned as follows:
(10)hl=[ξlηl],
where ξl∉I+l and ηl∈I+l. That is to say ξl is the active part of the current iteration. Then we compute the projection of ξ by the following equation:
(11)ξ=P[ξl-αD¯l∇g(ξl)],
where α≥0 is the iteration step-size, D¯l is a gradient scaling matrix, and P[a] is the positive projection, that is,
(12)P[a]={a,a>00,a≤0.
Next we update the current result hl by the following rules:
(13)hl+1=[ξηl]=[P[ξl-αD¯l∇g(ξl)]0],
where the right equation uses the fact that ηl is fixed to zeros, it can be comprehended as that this part satisfies the KKT optimal condition (8), then in the next iteration we need not to update this part. After updating the hl+1, we can compute the ∇g(hl+1) and update the inactive set I+l+1 to obtain the next iterative result. In the whole iteration we will adhere to the iteration rule until the result satisfies the stop criterions.


## 2.2. The Approximate toD¯l

As the size ofξl and ηl changes at each iteration, the computation of the matrix D¯l is not an easy task. Note that the curvature information of ξl is received from the curvature information of hl. Due to this fact, the gradient matrix D¯l can be obtained by taking the proper submatrix of the matrix D¯ to make this task easier, where D¯ is the matrix that carries the curvature information of the vector hl. In realization of the method we can try to eliminate the curvature information of ηl. Then the update rule of (13) can be regarded as follows:
(14)hl+1=P[hl-αD¯∇g(h)].
In [5], Kim et al. used the BFGS update method to approximate the Hessian matrix of the object function g(h). The BFGS method is well established and only uses the gradient information of the object function. But BFGS needs the time consumption. Many researchers have experimentally observed that the symmetric rank-one (SR1) rule performs better than BFGS quasi-Newton update rules [13]. In this paper, we will use the SR1 rule to approximate the Hessian matrix of g(h) to improve the the converge speed.Suppose thatHl is the current approximation of the Hessian matrix, then we apply the following SR1 rule to approximate the next Hessian matrix:
(15)Hl+1=Hl+(wl-Hlvl)(w-Hlvl)T(w-Hlvl)Twl,
where wl=∇g(hl+1)-∇g(hl) and vl=hl+1-hl. Let Dl be the inverse of the Hl. By using the Sherman-Morrison-Woodbury formula [14], we can obtain the inverse matrix Dl+1 of Hl+1 as follows:(16)Dl+1=Dl+(vl-Dlwl)(vl-Dlwl)T(vl-Dlwl)Twl.
Thus, we complete the approximate to D¯l in the iterative processes.


## 2.3. Line Search Strategy for Step-Size

For (13), α>0, we next consider how to find an α with the largest function reduction. For this purpose, we may choose the searching rule as follows:(17)α=argmin⁡α≥0g(P[hl-α∇g(hl)]).
To sum up Sections 2.1–2.3, we obtain the following symmetric rank-one NNLS algorithm (see Algorithm 1).Algorithm 1:The SR1 algorithm for NNLS.
Step  1. start data
Chooseh0≥0, and the maximum iterative number maxiter,
Setl=0;
Step  2.
Compute the inactive set by using the rule (9);
Step  3. (update the current solution)
Renew current solution by (13), where using (16) to approximate
the inverse of Hessian matrix and (17) to compute the stepsize;
Step  4.
ifthe current result reaches the stop criterion then stop, the current
solution is the final solution;
else go to Step  2 and Step  3.


## 3. A Symmetric Rank-One Quasi-Newton Method for NMF

In the former part, we have discussed the algorithm for NNLS and then in this part we are going to talk about the algorithm for NMF. Note that, in Section1, we have analyzed that the NMF problem can be resolved into several nonnegative least squares. Due to this property, we can apply the symmetric rank-one method to our NMF problem.


### 3.1. Applying the SR1 Directly to NNLS Subproblems

Firstly, for (4), we can do vectorization of the matrix H as follows:
(18)g(H)=12∥V-WH∥F2=12vec⁡(H)T(WTW⋱WTW)vec⁡(H)-tr⁡((WTV)TH)+12tr⁡(VTV),
where tr⁡(B) is the trace of matrix B. By using the quasi-Newton idea, the iteration of vec⁡(H) can be rewritten as
(19)vec⁡(Hl+1)=P[vec⁡(Hl)-αlD¯∇g(H)],
where αl is the step-size of the current iteration, D¯ is the gradient scaling matrix, and ∇g(H) is the gradient of (18).Equation (18) is actually a nonnegative least square problem; then we can apply the method discussed in the former part to this problem directly. First let the inactive set of (18) be defined as
(20)IHl={i∣0≤(vec⁡(H))i≤ϵl,∇g(H)i>0},
where ϵl can be obtained by the same method as the former part (see Section 2.1). In (19), the gradient scaling matrix Dl¯ is approximated by the symmetric rank-one quasi-Newton method. The step-size αl is done by the line search.For (3), the object function corresponding to W may be also rewritten as
(21)g(W)=12∥VT-HTWT∥F2=12vec⁡(H)T(HHT⋱HHT)vec⁡(WT)-tr⁡((HVT)TWT)+12tr⁡(VTV).
Thus, we can solve it by the same method as (18).


### 3.2. The Algorithm for NMF

Next, based on Section3.1 and the Algorithm 1, we give the symmetric rank-one quasi-Newton algorithm (SR1) for NMF as in Algorithm 2.Algorithm 2:The SR1 algorithm for NMF.
Initialization
Input matrixV, the maximum number of iteration: maxiter;
initial matrixW0;
the approximate rankk, 1≤k≤min⁡{m,n}.
fori=1: maxiter
(1)W=Wl,  Hold=Hl;
(2) Compute the matrixHnew by using the Algorithm 1, where the object function is (18);
(3)Hold=Hnew,  Hl+1=Hold;
(4)H=Hl+1,  Wold=Wl;
Compute the matrixWnew by using the Algorithm 1, where the object function is (21);
(5)Wold=Wnew,  Wl+1=Wold;
(6) if the stopping criteria is met, then break.
end forObviously, the computation of subproblems (3) and (4) is high cost. Each subproblem requires an iterative procedure, which can be regarded as a subiteration, and we need to compute the following gradients:(22)∇fH(W,H)=WTWH-WTV,∇fW(W,H)=WHHT-VHT.
In order to reduce the above computation, we can compute WTW, WTV and HTH, VHT in the outer iteration in our method.


### 3.3. Computational Complexity of Algorithms

According to numerical experiments in Section4, we find that our algorithm has a slight advantage compared with the faster algorithms such as the fast projection-based method (fnmae) in [5] and alternating nonnegative least squares using the projected Newton method (pnm) [11]. To demonstrate this problem more effectively, next, let us give a simple comparison on their computational complexity.(i)
For the algorithm SR1, it is very clear that it costsO(nmk) to compute the matrices WTW and WTV in the subproblem of NMF and it costs O(k2m2) to compute the gradient scaling matrix and gradient matrix, while the computation of step-size needs O(smk2), where s is the iteration number of the step-size. So the computational complexity of the whole method is as follows:
(23)O(nmk)+iter*(O(smk2)+O(k2m2)),
where the iter is the iteration number, that is, maxiter.(ii)
The computational complexity of pnm in [11] is
(24)O(nmp)+iter*(O(smk2)+O(m~s~)),
where s is the iteration number of the step-size and iter is the number of inner iteration. In this method the matrix D¯ is expressed by block diagonal matrix, and m~ is the number of the subblocks and the s~ is the average column number of the subblocks.(iii)
The computational complexity of fnmae in [5] is as follows:
(25)O(nmk)+iter*(O(smk2)+O(m3k3)),
where iter is the iteration number to solve the matrix H  and s is the iteration number of solving the step-size.


## 3.1. Applying the SR1 Directly to NNLS Subproblems

Firstly, for (4), we can do vectorization of the matrix H as follows:
(18)g(H)=12∥V-WH∥F2=12vec⁡(H)T(WTW⋱WTW)vec⁡(H)-tr⁡((WTV)TH)+12tr⁡(VTV),
where tr⁡(B) is the trace of matrix B. By using the quasi-Newton idea, the iteration of vec⁡(H) can be rewritten as
(19)vec⁡(Hl+1)=P[vec⁡(Hl)-αlD¯∇g(H)],
where αl is the step-size of the current iteration, D¯ is the gradient scaling matrix, and ∇g(H) is the gradient of (18).Equation (18) is actually a nonnegative least square problem; then we can apply the method discussed in the former part to this problem directly. First let the inactive set of (18) be defined as
(20)IHl={i∣0≤(vec⁡(H))i≤ϵl,∇g(H)i>0},
where ϵl can be obtained by the same method as the former part (see Section 2.1). In (19), the gradient scaling matrix Dl¯ is approximated by the symmetric rank-one quasi-Newton method. The step-size αl is done by the line search.For (3), the object function corresponding to W may be also rewritten as
(21)g(W)=12∥VT-HTWT∥F2=12vec⁡(H)T(HHT⋱HHT)vec⁡(WT)-tr⁡((HVT)TWT)+12tr⁡(VTV).
Thus, we can solve it by the same method as (18).


## 3.2. The Algorithm for NMF

Next, based on Section3.1 and the Algorithm 1, we give the symmetric rank-one quasi-Newton algorithm (SR1) for NMF as in Algorithm 2.Algorithm 2:The SR1 algorithm for NMF.
Initialization
Input matrixV, the maximum number of iteration: maxiter;
initial matrixW0;
the approximate rankk, 1≤k≤min⁡{m,n}.
fori=1: maxiter
(1)W=Wl,  Hold=Hl;
(2) Compute the matrixHnew by using the Algorithm 1, where the object function is (18);
(3)Hold=Hnew,  Hl+1=Hold;
(4)H=Hl+1,  Wold=Wl;
Compute the matrixWnew by using the Algorithm 1, where the object function is (21);
(5)Wold=Wnew,  Wl+1=Wold;
(6) if the stopping criteria is met, then break.
end forObviously, the computation of subproblems (3) and (4) is high cost. Each subproblem requires an iterative procedure, which can be regarded as a subiteration, and we need to compute the following gradients:(22)∇fH(W,H)=WTWH-WTV,∇fW(W,H)=WHHT-VHT.
In order to reduce the above computation, we can compute WTW, WTV and HTH, VHT in the outer iteration in our method.


## 3.3. Computational Complexity of Algorithms

According to numerical experiments in Section4, we find that our algorithm has a slight advantage compared with the faster algorithms such as the fast projection-based method (fnmae) in [5] and alternating nonnegative least squares using the projected Newton method (pnm) [11]. To demonstrate this problem more effectively, next, let us give a simple comparison on their computational complexity.(i)
For the algorithm SR1, it is very clear that it costsO(nmk) to compute the matrices WTW and WTV in the subproblem of NMF and it costs O(k2m2) to compute the gradient scaling matrix and gradient matrix, while the computation of step-size needs O(smk2), where s is the iteration number of the step-size. So the computational complexity of the whole method is as follows:
(23)O(nmk)+iter*(O(smk2)+O(k2m2)),
where the iter is the iteration number, that is, maxiter.(ii)
The computational complexity of pnm in [11] is
(24)O(nmp)+iter*(O(smk2)+O(m~s~)),
where s is the iteration number of the step-size and iter is the number of inner iteration. In this method the matrix D¯ is expressed by block diagonal matrix, and m~ is the number of the subblocks and the s~ is the average column number of the subblocks.(iii)
The computational complexity of fnmae in [5] is as follows:
(25)O(nmk)+iter*(O(smk2)+O(m3k3)),
where iter is the iteration number to solve the matrix H  and s is the iteration number of solving the step-size.


## 4. Numerical Experiments

In this part, we present some experiment results from both rand data and real world data and compare our method to the following six algorithms:(1)
fnmae: fast projection-based method in [5];(2)
pnm: alternating nonnegative least squares using the projected Newton method in [11];(3)
AS: the active set method in [9];(4)
alsq: the alternating nonnegative least squares method in [2];(5)
nnmalq: the alternating nonnegative least squares method in [15];(6)
nmf: the projected gradient method in [3];(7)
SR1: the symmetric rank-one quasi-Newton method in this paper.In our experiments, we initialize all the methods randomly and show plots of the relative error against the number of iteration or the time of iteration, where the relative error of approximation is∥V-WH∥F/∥V∥F.


### 4.1. Synthetic Data Experiments

By generating the synthetic data randomly, we test the random generated matrixV∈ℝn×m(n=200,m=40) with the approximate rank k=10 and another matrix V∈ℝn×m(n=2000,m=800) with the approximate rank k=5,10, and 20, respectively. In our experiments, we yield ten of these random matrices and then we take the average value of the ten experiment data as the final result.Figure1 shows the relative error of approximation against iteration number of random matrix V∈ℝ200×40 for the rank k=10 for different methods. As a whole the relative errors of approximation of the SR1, fnmae, and pnm are very similar after a few iterations, but after 80 iterations we can see a slight advantage of our method.Figure 1
Relative errors of approximation against iteration number for SR1, fnmae, pnm, alsq, and nnmalq with the matrixV∈ℝ200×40 for the approximate rank k=10.Figures2, 3, and 4 are the numerical results of the random matrix V∈ℝ2000×800 with the approximate rank k=5,10,and20, respectively. The figures in (a) show the comparisons among the SR1, fnmae, pnm, As, alsq, and nnmalq methods. To have a clear look at similar numerical results of SR1, pnm, and As methods, we present the comparisons of the three methods in (b). From these figures, we can learn that the relative error of approximation of SR1, fnmae, pnm, and As are also similar after a few iterations under the same condition. But our method SR1 shows more robustness than the other ones in all cases.Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=5.
(a)
(b)Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=10.
(a)
(b)Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=20.
(a)
(b)Figures5 and 6 are the plots of relative error against the running time. One can see that our method converges much faster than the nnmalq and alsq methods. Comparing with fnmae, our method converges as fast as fnmae at the beginning; after seconds our method shows a slight advantage in the left side of Figure 5.Relative errors of approximation against time for SR1 and fnmae with the matrixV∈ℝ200×40 in (a) and V∈ℝ2000×800 in (b) for the approximate rank k=10.
(a)
(b)Figure 6
Relative errors of approximation against time for SR1, fnmae, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=10.The above numerical results show that our symmetric rank-one quasi-Newton method improves the efficiency of the quasi-Newton method every iteration and cuts down the time of each iteration. Comparing with alsq and nnmalq, our SR1 method decreases faster in each iteration, while our method has a similar convergence with fnmae, As, and pnm. The nnmalq method decreases slowly in every iteration, so becomes less competitive. In addition, we find that the error of nmf method is much larger than the other above six methods mentioned in our experiment, so we did not plot its error curve in Figures1–6.


### 4.2. Application to Imaging Processing

NMF was originally motivated by Lee and Seung [1] in an image processing application. Many others have also considered NMF to image processing, face recognition application, model recognition application, signal processing application, and so forth. In this part, we are going to do some numerical experiments on image processing. Our experiment is done on four random chosen faces (http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html) and the size of each face image is 92 × 112 with 256 gray per pixel. We give out the reconstructed image of running 10, 20, and 50 iterations, respectively. In this experiment, we set the approximate rank k=14, and the images are taken randomly from five random initial matrices in each iteration. The first image in each row is the original image; others are the reconstructed images by SR1, fnmae, pnm, AS, alsq, nnmalq, and nmf methods, respectively.From Figures7 and 8, we can observe that, after 10 and 20 iterations, SR1 can reconstruct almost the same as fnmae, AS, and pnm procedures and much better than alsq, nnmalq, and nmf procedures. Figure 9 illustrates that the images obtained via SR1, As, fnmae and pnm methods have similar quality, and the images obtained by alsq are still vague. These pictures show that our method SR1 is also applied to the imaging processing and has a better reconstructing ability.Figure 7
Reconstructed images after 10 iterations.Figure 8
Reconstructed images after 20 iterations.Figure 9
Reconstructed images after 50 iterations.


### 4.3. Application to Text Clustering

In order to test its ability in the other cases, we are going to apply our method to text clustering problem and compare with the other methods. We show the numerical results of the above methods except for the nmf method on four text data sets, which are high dimensional and sparse. The text data is collected from newspapers  (http://www.shi-zhong.com/software/docdata.zip). In all the numerical experiments, we let the approximate rank k be the number of classes and we initialize the matrix randomly.We present the comparisons of the six methods mentioned above except for the nmf method in (a) of Figures10, 11, 12, and 13. (b) is the convergent specifics on the SR1, AS, and pnm methods. From Figures 10–13, we can known that, at the beginning several iterations our method is much better than pnm one, and decrease much faster than the other methods with the same random initial matrix. After five iterations, the convergence of these methods are more and more similar.Relative error of approximation against the number of iteration for the text data classic 300.
(a)
(b)Relative error of approximation against the number of iteration for the text data classic 3891.
(a)
(b)Relative error of approximation against the number of iteration for the text data hitech.
(a)
(b)Relative error of approximation against the number of iteration for the text data tr45.
(a)
(b)


## 4.1. Synthetic Data Experiments

By generating the synthetic data randomly, we test the random generated matrixV∈ℝn×m(n=200,m=40) with the approximate rank k=10 and another matrix V∈ℝn×m(n=2000,m=800) with the approximate rank k=5,10, and 20, respectively. In our experiments, we yield ten of these random matrices and then we take the average value of the ten experiment data as the final result.Figure1 shows the relative error of approximation against iteration number of random matrix V∈ℝ200×40 for the rank k=10 for different methods. As a whole the relative errors of approximation of the SR1, fnmae, and pnm are very similar after a few iterations, but after 80 iterations we can see a slight advantage of our method.Figure 1
Relative errors of approximation against iteration number for SR1, fnmae, pnm, alsq, and nnmalq with the matrixV∈ℝ200×40 for the approximate rank k=10.Figures2, 3, and 4 are the numerical results of the random matrix V∈ℝ2000×800 with the approximate rank k=5,10,and20, respectively. The figures in (a) show the comparisons among the SR1, fnmae, pnm, As, alsq, and nnmalq methods. To have a clear look at similar numerical results of SR1, pnm, and As methods, we present the comparisons of the three methods in (b). From these figures, we can learn that the relative error of approximation of SR1, fnmae, pnm, and As are also similar after a few iterations under the same condition. But our method SR1 shows more robustness than the other ones in all cases.Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=5.
(a)
(b)Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=10.
(a)
(b)Relative errors of approximation against iteration number for SR1, fnmae, pnm, AS, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=20.
(a)
(b)Figures5 and 6 are the plots of relative error against the running time. One can see that our method converges much faster than the nnmalq and alsq methods. Comparing with fnmae, our method converges as fast as fnmae at the beginning; after seconds our method shows a slight advantage in the left side of Figure 5.Relative errors of approximation against time for SR1 and fnmae with the matrixV∈ℝ200×40 in (a) and V∈ℝ2000×800 in (b) for the approximate rank k=10.
(a)
(b)Figure 6
Relative errors of approximation against time for SR1, fnmae, alsq, and nnmalq with the matrixV∈ℝ2000×800 for the approximate rank k=10.The above numerical results show that our symmetric rank-one quasi-Newton method improves the efficiency of the quasi-Newton method every iteration and cuts down the time of each iteration. Comparing with alsq and nnmalq, our SR1 method decreases faster in each iteration, while our method has a similar convergence with fnmae, As, and pnm. The nnmalq method decreases slowly in every iteration, so becomes less competitive. In addition, we find that the error of nmf method is much larger than the other above six methods mentioned in our experiment, so we did not plot its error curve in Figures1–6.


## 4.2. Application to Imaging Processing

NMF was originally motivated by Lee and Seung [1] in an image processing application. Many others have also considered NMF to image processing, face recognition application, model recognition application, signal processing application, and so forth. In this part, we are going to do some numerical experiments on image processing. Our experiment is done on four random chosen faces (http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html) and the size of each face image is 92 × 112 with 256 gray per pixel. We give out the reconstructed image of running 10, 20, and 50 iterations, respectively. In this experiment, we set the approximate rank k=14, and the images are taken randomly from five random initial matrices in each iteration. The first image in each row is the original image; others are the reconstructed images by SR1, fnmae, pnm, AS, alsq, nnmalq, and nmf methods, respectively.From Figures7 and 8, we can observe that, after 10 and 20 iterations, SR1 can reconstruct almost the same as fnmae, AS, and pnm procedures and much better than alsq, nnmalq, and nmf procedures. Figure 9 illustrates that the images obtained via SR1, As, fnmae and pnm methods have similar quality, and the images obtained by alsq are still vague. These pictures show that our method SR1 is also applied to the imaging processing and has a better reconstructing ability.Figure 7
Reconstructed images after 10 iterations.Figure 8
Reconstructed images after 20 iterations.Figure 9
Reconstructed images after 50 iterations.


## 4.3. Application to Text Clustering

In order to test its ability in the other cases, we are going to apply our method to text clustering problem and compare with the other methods. We show the numerical results of the above methods except for the nmf method on four text data sets, which are high dimensional and sparse. The text data is collected from newspapers  (http://www.shi-zhong.com/software/docdata.zip). In all the numerical experiments, we let the approximate rank k be the number of classes and we initialize the matrix randomly.We present the comparisons of the six methods mentioned above except for the nmf method in (a) of Figures10, 11, 12, and 13. (b) is the convergent specifics on the SR1, AS, and pnm methods. From Figures 10–13, we can known that, at the beginning several iterations our method is much better than pnm one, and decrease much faster than the other methods with the same random initial matrix. After five iterations, the convergence of these methods are more and more similar.Relative error of approximation against the number of iteration for the text data classic 300.
(a)
(b)Relative error of approximation against the number of iteration for the text data classic 3891.
(a)
(b)Relative error of approximation against the number of iteration for the text data hitech.
(a)
(b)Relative error of approximation against the number of iteration for the text data tr45.
(a)
(b)


## 5. Conclusion

In this paper, we present an algorithm for NMF problem, it is different from the famous fnmae method [5] in the following aspects.(1) The active set in our method is a relaxed form, while the active set in [5] is a hard form. The fact that an active set of hard form exhibits undesirable discontinuity at the boundary of the constraint set has been shown in Bertsekas [12], which is not good to the convergence rate, so we use the relax form to avoid this problem.(2) In addition, the fnmae method [5] uses the BFGS method to approximate the Hessian matrix. However, we apply the symmetric rank-one method to approximate the Hessian matrix. The symmetric rank-one method has been experimentally shown to be better than the BFGS method.What is more, our method is also different from the pnm method in [11]. The paper [11] uses the symmetric property of the Hessian matrix to approximate the inverse of the Hessian matrix by the Cholesky factorization, while in our paper we approximate the inverse of the Hessian matrix by the symmetric rank-one method. Moreover, some numerical experiments show that the object function decreases more per iteration in our method than in the pnm method [11].To summarize, we propose a symmetric rank-one quasi-Newton method on the NMF. This method maintains the decrease speed per iteration and decrease the computational time. From the experimental results both the synthetic data and real world data, our method performs very well, which shows more robustness than the other ones in almost all cases. In addition, there are many factors which influence the efficiency of the algorithm for NMF, such as the approximate rank, the stop criterion, and sparseness. On these problems, we will continue researching in the future. Whats more we will also research its usage in the pattern recognition in our next research.



---

*Source: 846483-2014-01-22.xml*