# Two Efficient Algorithms for Orthogonal Nonnegative Matrix Factorization

Authors: Jing Wu; Bin Chen; Tao Han
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Nonnegative matrix factorization (NMF) is a popular method for the multivariate analysis of nonnegative data. It involves decomposing a data matrix into a product of two factor matrices with all entries restricted to being nonnegative. Orthogonal nonnegative matrix factorization (ONMF) has been introduced recently. This method has demonstrated remarkable performance in clustering tasks, such as gene expression classification. In this study, we introduce two convergence methods for solving ONMF. First, we design a convergent orthogonal algorithm based on the Lagrange multiplier method. Second, we propose an approach that is based on the alternating direction method. Finally, we demonstrate that the two proposed approaches tend to deliver higher-quality solutions and perform better in clustering tasks compared with a state-of-the-art ONMF.


---

## Body

## 1. Introduction

Nonnegative matrix factorization (NMF) has been investigated by many researchers, such as Paatero and Tapper [1]. However, it was only through the works of Lee and Seung published in Nature and NIPS [2, 3] that this method has gained popularity. Based on the argument that nonnegativity is crucial for human perception, they proposed a simple algorithm to find nonnegative representations of nonnegative data and images. The basic NMF is a technique that decomposes a nonnegative data matrix into a pair of other nonnegative matrices with lower ranks: (1)Y≈AX,where Y∈ℝ+I×T denotes a nonnegative data matrix. A∈ℝ+I×J and X∈ℝ+J×T denote the basis matrix and coefficient matrix, respectively. J denotes the number of factors that is usually chosen so that J≪minI,T. Recently, NMF has proven to be useful for many applications in pattern recognition, multimedia, text mining, and clustering, as well as environment and DNA gene expression classifications [4–12].In solving the NMF problem, the update rules given by Lee and Seung [3] take a multiplicative form and satisfy the nonnegative constraint implicitly and elegantly. The extension of the multiplicative updates provided by them [3] can be found, for instance, in the works of Sha et al. [10], whereby a multiplicative update rule is proposed to solve a nonnegative quadratic programming problem. However, the slow convergence of the multiplicative updating rule has been pointed out, and more efficient algorithms equipped with more powerful theoretical convergence properties have been introduced. These efficient algorithms are based on either the alternating nonnegative least squares (ANLS) framework or the hierarchical alternating least squares (HALS) method.Recently, other forms of update rules have been studied to increase the speed of convergence. For instance, a quasi-Newton method for NMF was considered by Zdunek and Cichocki [13], and a projected gradient method was proposed by Hoyer and Lin [6, 14].There has also been significant interest in orthogonal nonnegative matrix factorization (ONMF), which imposes orthogonality constraints on the factor matrix [15–18]. It has been proven that ONMF is equivalent to K-means clustering [16, 18]. Moreover, as a result of the nonnegativity and orthogonality constraints, the orthogonal factor matrix is naturally sparse. In such cases, the factorization is essentially unique, and we can ignore the permutation problem.Ding et al. and Pompili et al. [16, 18] proposed two different kinds of methods to solve the ONMF problem. In a study conducted by Ding et al. [16], ONMF algorithms strictly enforce nonnegativity for each iterate while trying to achieve orthogonality at the limit. This can be done using a proper penalization term [16], a projection matrix formulation [19], or choosing a suitable search direction [15]. In a study conducted by Pompili et al. [18], the authors proposed a method working the opposite way: at each iteration, a projected gradient scheme is used to ensure that the orthogonal factor iterates are orthogonal but not necessarily nonnegative. However, these algorithms still are not convergent algorithms for ONMF. In this study, we first briefly introduce NMF and ONMF algorithms. Next, based on the Lagrange multiplier, we design a convergence algorithm, which is a method similar to the algorithm presented by Lin [20]. Because NMF algorithms based on the alternating direction method are more efficient than multiplicative update algorithms, we propose another convergence algorithm for solving the ONMF problem by combining the ADM approach with our convergence algorithm. Experiments on two grayscale images and the ALLAML gene-sample data demonstrate that Algorithms 2 and 3 tend to deliver higher -quality solutions and perform better in a clustering task compared with Algorithm 1.Algorithm 1: AUR algorithm for the ONMF.
Input:Y∈ℝ+I×T and J∈ℕ.Output: FactorsA∈ℝ+I×J,X∈ℝ+J×T, which solve problem (18).(1)
Select initial matrixA0>0 and X0>0.(2)
Fork = 1 : maxIt(2.1)
Calculate the factor matrixA using the following formula:aij←aij−a¯ij/AXXTij+δ∇ADFAij(2.2)
Calculate the factor matrixX using the following formula:xjt←xjt−x¯jt/ATAX+αXXTXjt+δ∇XDFXjtEnd DoAlgorithm 2: Modified AUR algorithm of ONMF.
Input:Y∈ℝ+I×T and J∈ℕ.Output:A∈ℝ+I×J,X∈ℝ+J×T, which solve the optimization problem (1.1).(1)
Select initial matricesA0>0 and X0>0.(2)
Fork = 1:maxIt(2.1)
Calculate the factor matrixA using the following formula:aij←maxaij−a¯ij/AXXTij+δ∇ADFAij,ε(2.2)
Calculate the factor matrixX using the following formula:xjt←maxxjt−x¯jt/ATAX+αXXTX+ΛYjt+δ∇XDFXjt,ε(2.3)
Calculate the factor matrixΛ using the following formula:Λ←Λ+γρXXT−IEnd DoAlgorithm 3: ADM-based algorithm of ONMF.
Input:Y∈ℝ+I×T and J∈ℕ.Output:A∈ℝ+I×J,X∈ℝ+J×T, which solve the optimization problem (1.1).(1)
Select initial matricesA0>0 and X0>0.(2)
Fork = 1:maxItUpdateAk,Xk,Uk,Vk,Λk,Πk by formulas (62)–(67)If a stopping criterion is met, thenexit and outputAk,Xkend ifEnd DoThe remainder of this article is organized as follows: Section2 gives a brief review of NMF and ONMF algorithms. In Section 3, an improved version of the original ONMF algorithm and the convergence analysis is provided and another new ONMF algorithm based on the alternating least squares (ALS) approach is provided. Section 4 is devoted to numerical experiments. Some concluding remarks are provided in Section 5.


## 2. Related Algorithms

In this section, we provide a brief review of NMF and ONMF algorithms.


### 2.1. Nonnegative Matrix Factorization

NMF aims to find a nonnegative matrixY=y1,y2,…,yr∈ℝ+m×r and another nonnegative matrix XT=B=b1,b2,…,br∈ℝ+n×r whose product approximate a given nonnegative matrix Y∈ℝ+m×n:(2)minA≥0,X≥0Y−ABTF2.Because the NMF problem is not convex in bothA and B, various algorithms have been proposed. The first method for solving the equation (2) is the ALS algorithm utilized by Paatero and Tapper in 1994. It minimizes the least squares cost function with respect to either A or B, one at a time, by fixing the other and disregarding nonnegativity; after each least squares step, the algorithm sets all negative entries to zero. The scheme can be updated as follows:(3)A←YXTXXT†,(4)X←ATA†ATY,where A† denotes the Moore–Penrose pseudoinverse.Another popular method for NMF is the multiplicative updating method proposed by Lee and Seung. It can be expressed as follows:(5)A←A⊙YXT⊘AXXT,(6)X←X⊙ATY⊘ATAX,where ⊙ and ⊘ denote elementwise multiplication and division, respectively. When they are started from nonnegative initial estimates, iterates will remain nonnegative throughout the iterations. By all indications, this algorithm appears to be the most widely used solution method in NMF by far.The modified update rule given in equations (7) and (8) was first proposed by Gillis and Glineur [21]; it can be expressed as follows:(7)A←maxɛ,A⊙YXT⊘AXXT,(8)X←maxɛ,X⊙ATY⊘ATAX,where ɛ>0 denotes a small number. They proved that if a sequence of solutions generated by (7) and (8) has a limit point, then this point is necessarily a stationary point of the following optimization problem:(9)minA≥ε,X≥εY−ABTF2.


### 2.2. Original ONMF Algorithms

In an ONMF, an additional orthogonal constraint,XXT=I, is imposed on the NMF. We briefly review the first ONMF algorithm proposed by Ding et al. [16] and reveal the problem behind ONMF.The goal of ONMF is to find a nonnegative orthogonal matrixX and a nonnegative matrix A by minimizing the following objective function using a Lagrangian multiplier λ. Because it is not easy to solve this problem for every value of λ, Ding et al. [16] ignored the nonnegativity and relied on XXT=I to obtain a unique value of ΛX:(10)LA,X=12Y−AXF2+12trΛXXXT−I.The KKT conditions for the objective in equation (10) are as follows:(11)A∗≥0,X∗≥0,∇AA∗≥0,∇XX∗≥0,∇AA∗⊙A∗=0,∇XX∗⊙X∗=0,with(12)∇AA=AXXT−YXT,∇XX=ATAX−ATY+ΛXX.Using the multiplicative updating method, the update rules ofA and X are derived as follows:(13)aij←aijYXTijAXXTij,xjt←xjtATYjtATA+ΛXXjt.The problem with this algorithm is that it is difficult to determineΛc. By summing over the index i, the authors found an exact formulation for the diagonal entries:(14)Λcii=ATYXT−ATAii.The off-diagonal entries are obtained by ignoring the nonnegativity constraint onX and setting ∇XX to a zero matrix:(15)∇XX=ATAX−ATY+ΛXX=0.(16)ΛXij=ATYXT−ATAij,∀i≠j.Equation (15) is derived from equation (10) using the fact that XF2=trXTX and from equation (16) using the orthogonality constraint XXT=I. By combining equations (14) and (16), ΛX can be defined as follows:(17)Λc=ATYXT−ATA.Note that the formulation with the specific value ofΛX does not strictly satisfy orthogonality. Therefore, ONMF algorithms prioritize the approximation while relaxing the degree of orthogonality.To improve the approximation, Mirzal [22] proposed the following cost function:(18)minA,XDFYAX=12Y−AXF2+α2XXT−IF2.Using the same strategy as in the aforementioned algorithm, the update forms are derived as follows:(19)aij←aijYXTijAXXTij,xjt←xjtATY+αXjtATAX+αXXTXjt.Next, the author proposed an equivalent algorithm with a robust convergence guaranteed by (1) transforming the algorithm into a corresponding AUR algorithm and (2) replacing every zero entry that does not satisfy the KKT conditions with a small positive number to prevent zero locking. This strategy was employed to derive a convergent algorithm for the ONMF.The AUR version of equation (18) can be expressed as follows:(20)aij←aij−aijAXXTij∇ADFAij,xjt←xjt−xjtATAX+αXXTXjt∇XDFXjt,where Aij=aij denotes the ijth entry of matrix A. Additionally, the author proposed an equally robust AUR algorithm (see Algorithm 1) for equation (18), where maxIt denotes a limit maxtime of iteration, ∇ADF represents the gradient of the objective function (18), and(21)a¯ijk=aijk,if∇AAk,Xkij≥0,maxaijk,δ,if∇AAk,Xkij<0,x¯jtk=xjtk,if∇XAk+1,Xkjt≥0,maxxijk,δ,if∇XAk+1,Xkjt<0.The author proved the convergence of this algorithm. However, this algorithm requires a high computational cost. The algorithm that is based on the alternated application of these rules in equations (5) and (6) is not guaranteed to converge to a first-order stationary point. However, a slight modification proposed by Lin [20] achieves this property (generally speaking, MU is recast as a variable metric steepest descent method, and the step length is modified accordingly). In another study [23], the authors demonstrated through numerical experiments that the update rules in equations (7) and (8) work better than the original update rules in equations (5) and (6) in some cases. This indicates that equations (7) and (8) are important not only from a theoretical point of view but also from a practical viewpoint. Therefore, in the next section, we propose convergence updates based on the alternated approach.


## 2.1. Nonnegative Matrix Factorization

NMF aims to find a nonnegative matrixY=y1,y2,…,yr∈ℝ+m×r and another nonnegative matrix XT=B=b1,b2,…,br∈ℝ+n×r whose product approximate a given nonnegative matrix Y∈ℝ+m×n:(2)minA≥0,X≥0Y−ABTF2.Because the NMF problem is not convex in bothA and B, various algorithms have been proposed. The first method for solving the equation (2) is the ALS algorithm utilized by Paatero and Tapper in 1994. It minimizes the least squares cost function with respect to either A or B, one at a time, by fixing the other and disregarding nonnegativity; after each least squares step, the algorithm sets all negative entries to zero. The scheme can be updated as follows:(3)A←YXTXXT†,(4)X←ATA†ATY,where A† denotes the Moore–Penrose pseudoinverse.Another popular method for NMF is the multiplicative updating method proposed by Lee and Seung. It can be expressed as follows:(5)A←A⊙YXT⊘AXXT,(6)X←X⊙ATY⊘ATAX,where ⊙ and ⊘ denote elementwise multiplication and division, respectively. When they are started from nonnegative initial estimates, iterates will remain nonnegative throughout the iterations. By all indications, this algorithm appears to be the most widely used solution method in NMF by far.The modified update rule given in equations (7) and (8) was first proposed by Gillis and Glineur [21]; it can be expressed as follows:(7)A←maxɛ,A⊙YXT⊘AXXT,(8)X←maxɛ,X⊙ATY⊘ATAX,where ɛ>0 denotes a small number. They proved that if a sequence of solutions generated by (7) and (8) has a limit point, then this point is necessarily a stationary point of the following optimization problem:(9)minA≥ε,X≥εY−ABTF2.


## 2.2. Original ONMF Algorithms

In an ONMF, an additional orthogonal constraint,XXT=I, is imposed on the NMF. We briefly review the first ONMF algorithm proposed by Ding et al. [16] and reveal the problem behind ONMF.The goal of ONMF is to find a nonnegative orthogonal matrixX and a nonnegative matrix A by minimizing the following objective function using a Lagrangian multiplier λ. Because it is not easy to solve this problem for every value of λ, Ding et al. [16] ignored the nonnegativity and relied on XXT=I to obtain a unique value of ΛX:(10)LA,X=12Y−AXF2+12trΛXXXT−I.The KKT conditions for the objective in equation (10) are as follows:(11)A∗≥0,X∗≥0,∇AA∗≥0,∇XX∗≥0,∇AA∗⊙A∗=0,∇XX∗⊙X∗=0,with(12)∇AA=AXXT−YXT,∇XX=ATAX−ATY+ΛXX.Using the multiplicative updating method, the update rules ofA and X are derived as follows:(13)aij←aijYXTijAXXTij,xjt←xjtATYjtATA+ΛXXjt.The problem with this algorithm is that it is difficult to determineΛc. By summing over the index i, the authors found an exact formulation for the diagonal entries:(14)Λcii=ATYXT−ATAii.The off-diagonal entries are obtained by ignoring the nonnegativity constraint onX and setting ∇XX to a zero matrix:(15)∇XX=ATAX−ATY+ΛXX=0.(16)ΛXij=ATYXT−ATAij,∀i≠j.Equation (15) is derived from equation (10) using the fact that XF2=trXTX and from equation (16) using the orthogonality constraint XXT=I. By combining equations (14) and (16), ΛX can be defined as follows:(17)Λc=ATYXT−ATA.Note that the formulation with the specific value ofΛX does not strictly satisfy orthogonality. Therefore, ONMF algorithms prioritize the approximation while relaxing the degree of orthogonality.To improve the approximation, Mirzal [22] proposed the following cost function:(18)minA,XDFYAX=12Y−AXF2+α2XXT−IF2.Using the same strategy as in the aforementioned algorithm, the update forms are derived as follows:(19)aij←aijYXTijAXXTij,xjt←xjtATY+αXjtATAX+αXXTXjt.Next, the author proposed an equivalent algorithm with a robust convergence guaranteed by (1) transforming the algorithm into a corresponding AUR algorithm and (2) replacing every zero entry that does not satisfy the KKT conditions with a small positive number to prevent zero locking. This strategy was employed to derive a convergent algorithm for the ONMF.The AUR version of equation (18) can be expressed as follows:(20)aij←aij−aijAXXTij∇ADFAij,xjt←xjt−xjtATAX+αXXTXjt∇XDFXjt,where Aij=aij denotes the ijth entry of matrix A. Additionally, the author proposed an equally robust AUR algorithm (see Algorithm 1) for equation (18), where maxIt denotes a limit maxtime of iteration, ∇ADF represents the gradient of the objective function (18), and(21)a¯ijk=aijk,if∇AAk,Xkij≥0,maxaijk,δ,if∇AAk,Xkij<0,x¯jtk=xjtk,if∇XAk+1,Xkjt≥0,maxxijk,δ,if∇XAk+1,Xkjt<0.The author proved the convergence of this algorithm. However, this algorithm requires a high computational cost. The algorithm that is based on the alternated application of these rules in equations (5) and (6) is not guaranteed to converge to a first-order stationary point. However, a slight modification proposed by Lin [20] achieves this property (generally speaking, MU is recast as a variable metric steepest descent method, and the step length is modified accordingly). In another study [23], the authors demonstrated through numerical experiments that the update rules in equations (7) and (8) work better than the original update rules in equations (5) and (6) in some cases. This indicates that equations (7) and (8) are important not only from a theoretical point of view but also from a practical viewpoint. Therefore, in the next section, we propose convergence updates based on the alternated approach.


## 3. Alternating Direction Algorithm for ONMF

### 3.1. Classic ADM Approach

ADM is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the multipliers method. The algorithm solves problems in the form.(22)minfx+gzsubject toAx+Bz=c,with variables x∈ℝn and z∈ℝm, where A∈ℝp×n, B∈ℝp×m, and c∈ℝp. We assume that f and g are convex. The optimal value of equation (22) will be denoted by(23)p∗=inffx+gz|Ax+Bz=c.As in the method of ADM, we form the augmented Lagrangian as follows:(24)Lρx,z,y=fx+gz+yTAx+Bz−c+ρ2Ax+Bz−c22.ADM consists of the following iterations:(25)xk+1≔argminxLρx,zk,yk,(26)zk+1≔argminzLρxk+1,z,yk,(27)yk+1≔yk+γρAxk+1+Bzk+1−c,where ρ>0 and γ∈0,1.618 represent a step length. The algorithm is similar to dual ascent and the method of multipliers. It consists of an x-minimization step (26), a z-minimization step (27), and a dual variable update (27). As in the method of multipliers, the dual variable update uses a step size that is equal to the augmented Lagrangian parameter ρ. Conversely, the classic augmented Lagrangian multiplier method requires a joint minimization with respect to both x and z, i.e., it replaces steps (26) and (27) with(28)xk+1,zk+1←argminx,zLAx,z,yT,which involves both fx and gz and could become much more expensive. In Section 3.2, we propose improved algorithms for ONMF and present the convergence analysis. We discuss the orthogonality constraint on the rows of X. Similar results can be derived for A.


### 3.2. Classic ADM Approach Extension to ONMF

First, we define the objective of the ONMF as follows:(29)minA,XLA,X=12Y−AXF2+α2XXT−IF2+Λ⋅XXT−I,such that A≥0andX≥0.The AUR version of Problem (18) can be modified as follows:(30)aij←aij−a¯ijAXXTij+δ∇ADFAij,xjt←xjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt.Inspired by the approach proposed by Gillis and Glineur [23] and Yoo and Choi [24], we employ these strategies to devise a convergent algorithm for NMF with an orthogonality constraint, which is Problem (4.1). The modified AUR version of Problem (4.1) can be modified as follows:(31)aij←maxaij−a¯ijAXXTij+δ∇ADFAij,ε,xjt←maxxjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt,ε,Λ←Λ+γρXXT−I,where ε>0 denotes a small number and ρ>0 and γ∈0,1.618 represent a step length.(32)a¯ijk=aijk,if∇AAk,Xkij≥0,maxaijk,σ,otherwise,x¯jtk=xjtk,if∇XAk+1,Xkjt≥0,maxxijk,σ,otherwise,stand for the modifications for avoiding zero locking with a small positive number σ and(33)∇AA=AXXT−YXT,∇XX=ATAX−ATY+αXXTX−αX+ΛY.The detailed procedure for the ALS-based version of Problem (4.1) can be seen in Algorithm2.For Algorithm2, we obtain the following basic result:Lemma 1.
For any natural numberk, if A0≥0 and X0≥0, then Ak>0 and Xk>0.Proof.
This statement forA is similar to the proof for X. Therefore, we only prove the statement for X. It is obvious for k=0; hence, we only need to prove the results for k≥1.Case 1.
For∇XAk+1,Xkjt≥0,(34)xjtk+1=ATAXk+αXkXkTXk+ΛYjtxjtk+δXkxjtkATAXk+αXkXkTXk+ΛYjt+δXk−ATAXk+αXkXkTXk+ΛYjtxjtkATAXk+αXkXkTXk+ΛYjt+δXk+Ak+1TY+αXkjtxjtkATAXk+αXkXkTXk+ΛYjt+δXk=δXk+Ak+1TY+αXkjtxjtkAk+1TYXkTjt+δXk.
With the nonnegativity ofAk+1T,X, and Y,(35)xjtk≥0⇒xjtk+1≥0.
Because(36)xjt←maxxjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt,ε,(37)xjtk+1>0.Case 2.
For∇XAk+1,Xkjt<0,(38)xjtk+1=xjtk−maxxijk,δ∇XAk+1,Xkjt∇ADFAij++δXk.
Becausemaxxijk,δ>0 and ∇XAk+1,Xkjt<0, equation (27) holds.
Similarly, the statement ofA can be proven.


#### 3.2.1. Convergence Analysis

To analyze the convergence of Algorithm2, we first prove its nonincreasing property under the update rules in equations (25) and (26), i. e.,(39)LAk+1,Xk+1≤LAk,Xk.Because the algorithm solves Problem (1.1) in an alternating method,LAk and LXk can be analyzed separately. Because the update rule of A is similar to that of X, it is sufficient to prove the nonincreasing property of LXk.By making use of the auxiliary function approach [3], the nonincreasing property of LXk can be proved by showing that(40)LXk+1≤GAk+1,Xk+1≤GAk+1,Xk≤GAk,Xk≤LXk,where the auxiliary function G is defined in a similar function by [22](41)GX,Xk=LXk+trX−XkT∇XLXk+12trX−XkTDX−Xk,(42)X=x1x2⋱xT∈ℝ+TJ×J,where xii=1,2,…,T denotes the ith column of X. Meanwhile, we define(43)∇XLXk=∇XLXk1∇XLXk2⋱∇XLXkT∈ℝ+TJ×J,where ∇XLXki denotes the ith column of ∇XLXk, and(44)D=diagD1,D2,…,DT∈ℝ+TJ×TJ,where Dtt=1,…,T denotes a diagonal matrix with its diagonal entries(45)djjt=ATAX+αXXTX+ΛYjt+δXkx¯jtk,ifj∈ln,0,otherwise.Here,ln=j|xjtk>0,∇XLXkjt≠0,orxjtk=0,∇XLXkjt<0 denotes the set of non-KKT indices in the nth column of Xk.Using the Taylor series, the expansion formulation forLX can be expressed as follows:(46)LX=LXk+trX−XkT∇XLXk+12trX−XkTLX−Xk+ξXk,where L=∇X2LXk and ξXk represent the higher components of the Taylor series, and(47)∇X2LXk=∇X2LXk1∇X2LXk2⋱∇X2LXkT∈ℝ+TJ×TJ,where ∇X2LXkn denotes the nth column of ∇X2LXk. Next, for the function GX,Xk, which is an auxiliary of LX, we first have to prove the following Lemma 2. As shown later, Ak and Xk must be bounded for LAk,Xk. The boundedness of Ak and Xk will be proven in Theorem 3.Lemma 2.
Given sufficiently largeδXk and the boundedness of Ak and, then(1)
GX,X=LX.(2)
GXk,Xk=LXk.(3)
GX,X≤GX,Xk. Moreover, GX,X=GX,Xk only if Xk satisfies the KKT conditions.(4)
GX,Xk≤GXk,Xk, and the equality holds only if Xk satisfies the KKT conditions.Proof.
(1) and (2) are obvious from equation (41).
From (3), becauseGX,X=LX, we have(48)GX,Xk−GX,X=GX,X−LX=12trX−XkTD−LX−Xk−ξXk=∑t=1Txt−xtkTDt−Ltx−xtk−ξXk,(49)Dt−Lt=D¯t+δXkD^t−Lt,where D¯t and δXkD^t are diagonal matrices that sum up to Dt with(50)d¯jjt=∇XDFXjt+x¯jtk,ifj∈ln,0,otherwise,d^jjt=1x¯jtk,ifj∈ln,0,otherwise.
Letxt−xtk=x^t; then,(51)GX,Xk−LX=12∑t=1T∑j=1Jx^jt2d¯jjt+δXk∑j=1Jx^jt2d^jjt−12∑t=1Tx^tkTLtx^t−ξXk.
With the boundedness ofAk and Xk and using sufficiently large δXk, the inequality GX,X≤GX,Xk can be guaranteed.
Additionally, it is clear thatGX,X=GX,Xk only if Xk satisfies the KKT conditions demonstrated in equation (48) regardless of ξXk. Because δXk is a variable, GX,X=GX,Xk holds only if X=Xk in equation (51) as a result of the update rule of X and the boundedness of Ak and Xk.
(4) Using(52)GXk,Xk−GX,Xk=−trX−XkT∇XLXkX−Xk−12trX−XkTDX−Xk,and the positive semidefinite of matrix D, GX,Xk is a strictly convex function with respect to Xk; it has a unique minimum that satisfies(53)DX−Xk+∇XLXk=0.
The aforementioned consideration results in the following:(54)GXk,Xk−GX,Xk=12trX−XkTDX−Xk≥0.
IfXk satisfies the KKT conditions, then X=Xk, and therefore, GXk,Xk=GX,Xk.
Conversely, assuming that the equality holds, butXk does not satisfy the KKT conditions, there is at least an index such that(55)xjt≠xjtk,∇XDFXjt++δXkx¯jtk≥δXkx¯jtk.
Furthermore, ifx¯jtk=0, then xjt=xjtk, which contradicts the equality assumption. Therefore, x¯jtk≠0. As a result,(56)GXk,Xk−GX,Xk≥xjt−xjtkδXkx¯jtk>0,which contradicts the assumption.Theorem 1.
Given a sufficiently large numberδXk and supposing that Ak and Xk are bounded, then for the auxiliary function G,(1)
LXk+1≤LXk under the update rule in equation (4) and the equality holds only if Xk satisfies the KKT conditions(2)
LAk+1≤LAk under the update rule in equation (3) and the equality holds only if Ak satisfies the KKT conditions(3)
LAk+1,Xk+1≤LAk,Xk and the equality holds only if Ak,Xk is a stationary pointProof.
(1) and (2) can be obtained directly from Lemma1. (3) By combining (1) and (2), it is obvious that(57)LAk+1,Xk+1≤LAk,Xk,with LAk+1,Xk+1=LAk,Xk only if Ak and Xk satisfy the KKT conditions, i.e., Ak,Xk is a stationary point.Theorem 2.
The objective function (22) is nonincreasing under the multiplicative update rules in equations (25) and (26).Proof.
This theorem is the corollary of Theorem1.Lemma 3.
SequenceAk,Xk,Λk generated by Algorithm 2 is bounded and LAk,Xk is lower bounded.Proof.
Let us verify Assumptions A1–A3 (assumptions A1–A3 can be found in the works of [23]). Assumption A1 holds because of the coercivity of ⋅F. Assumptions A2 and A3 hold because all the coefficient matrices are identity matrices. Additionally, because ⋅F is Lipschitz differentiable, Ak,Xk,Λk generated by Algorithm 3 is bounded, and LAk,Xk is lower bounded according to Theorem 2 presented in the works of [23].
The convergence property of the modified update rules in equations (27) and (29) is stated as follows:Theorem 3.
Given a sufficiently large numberδXk, every limit point of this algorithm is a stationary point of the following optimization problem:(58)minA≥ε,X≥εLA,X=12Y−AXF2+α2XXT−IF2+Λ⋅XXT−I.Proof.
BecauseAk,Xkk=1,…,∞ are in a compact set, the sequence Ak,Xk generated by Algorithm 2 has at least one limit point. Gillis and Glineur [23] introduced a similar analysis that every limit point of this algorithm is a stationary point of the optimization problem (41).


### 3.3. Auxiliary Variable-Based ADM Extension to ONMF

To facilitate the efficient use of ADM, we first introduce two auxiliary variables,U and V, and consider the following equivalent model:(59)minX,Y,U,V12Y−AXF2s.t.A=U,X=V,U≥0,V≥0,XXT=I,where A,U∈ℝm×q and Y,V∈ℝq×n. The augmented Lagrangian of equation (59) is as follows:(60)LAA,X,U,V,Λ,Π=12Y−AXF2+α2XXT−IF2+Λ⋅A−U+Π⋅X−V+ρ12A−UF2+ρ22X−VF2,where Λ∈ℝm×q,Π∈ℝq×n are Lagrangian multipliers, α,ρ1,ρ2>0 are penalty parameters, and A⋅B≔∑i,jaijbij for matrices A and B of the same size.The alternating direction method for equation (59) is derived by successively minimizing LA with respect to A,X,U,V, one at a time, while fixing others at their more recent values, i. e.,(61)Ak+1=argminLAA,Xk,Uk,Vk,Λk,Πk,Xk+1=argminLAAk+1,X,Uk,Vk,Λk,Πk,Uk+1=argminU≥0LAAk+1,Xk+1,U,Vk,Λk,Πk,Vk+1=argminU≥0LAAk+1,Xk+1,Uk+1,V,Λk,Πk,and then updating the multipliers Λ and Π. The introduction of the two auxiliary variables, U and V, makes it easy to carry out each of the alternating minimization steps. Specifically, these steps can be expressed in closed form as follows:(62)aij←aijYXT+ρ1UijAXXT+Π+ρ1Aij+δω,(63)xjt←xjtATY+αX+ρ2VjtATAX+αXXTX+Λ+ρ2Xjt+δω,(64)U←P+A+Λρ1,(65)V←P+X+Πρ2,(66)Λ←Λ+γρ1A−U,(67)Π←Π+γρ2X−V.The detailed procedure for the algorithm of Problem59 can be seen in Algorithm 3.


## 3.1. Classic ADM Approach

ADM is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the multipliers method. The algorithm solves problems in the form.(22)minfx+gzsubject toAx+Bz=c,with variables x∈ℝn and z∈ℝm, where A∈ℝp×n, B∈ℝp×m, and c∈ℝp. We assume that f and g are convex. The optimal value of equation (22) will be denoted by(23)p∗=inffx+gz|Ax+Bz=c.As in the method of ADM, we form the augmented Lagrangian as follows:(24)Lρx,z,y=fx+gz+yTAx+Bz−c+ρ2Ax+Bz−c22.ADM consists of the following iterations:(25)xk+1≔argminxLρx,zk,yk,(26)zk+1≔argminzLρxk+1,z,yk,(27)yk+1≔yk+γρAxk+1+Bzk+1−c,where ρ>0 and γ∈0,1.618 represent a step length. The algorithm is similar to dual ascent and the method of multipliers. It consists of an x-minimization step (26), a z-minimization step (27), and a dual variable update (27). As in the method of multipliers, the dual variable update uses a step size that is equal to the augmented Lagrangian parameter ρ. Conversely, the classic augmented Lagrangian multiplier method requires a joint minimization with respect to both x and z, i.e., it replaces steps (26) and (27) with(28)xk+1,zk+1←argminx,zLAx,z,yT,which involves both fx and gz and could become much more expensive. In Section 3.2, we propose improved algorithms for ONMF and present the convergence analysis. We discuss the orthogonality constraint on the rows of X. Similar results can be derived for A.


## 3.2. Classic ADM Approach Extension to ONMF

First, we define the objective of the ONMF as follows:(29)minA,XLA,X=12Y−AXF2+α2XXT−IF2+Λ⋅XXT−I,such that A≥0andX≥0.The AUR version of Problem (18) can be modified as follows:(30)aij←aij−a¯ijAXXTij+δ∇ADFAij,xjt←xjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt.Inspired by the approach proposed by Gillis and Glineur [23] and Yoo and Choi [24], we employ these strategies to devise a convergent algorithm for NMF with an orthogonality constraint, which is Problem (4.1). The modified AUR version of Problem (4.1) can be modified as follows:(31)aij←maxaij−a¯ijAXXTij+δ∇ADFAij,ε,xjt←maxxjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt,ε,Λ←Λ+γρXXT−I,where ε>0 denotes a small number and ρ>0 and γ∈0,1.618 represent a step length.(32)a¯ijk=aijk,if∇AAk,Xkij≥0,maxaijk,σ,otherwise,x¯jtk=xjtk,if∇XAk+1,Xkjt≥0,maxxijk,σ,otherwise,stand for the modifications for avoiding zero locking with a small positive number σ and(33)∇AA=AXXT−YXT,∇XX=ATAX−ATY+αXXTX−αX+ΛY.The detailed procedure for the ALS-based version of Problem (4.1) can be seen in Algorithm2.For Algorithm2, we obtain the following basic result:Lemma 1.
For any natural numberk, if A0≥0 and X0≥0, then Ak>0 and Xk>0.Proof.
This statement forA is similar to the proof for X. Therefore, we only prove the statement for X. It is obvious for k=0; hence, we only need to prove the results for k≥1.Case 1.
For∇XAk+1,Xkjt≥0,(34)xjtk+1=ATAXk+αXkXkTXk+ΛYjtxjtk+δXkxjtkATAXk+αXkXkTXk+ΛYjt+δXk−ATAXk+αXkXkTXk+ΛYjtxjtkATAXk+αXkXkTXk+ΛYjt+δXk+Ak+1TY+αXkjtxjtkATAXk+αXkXkTXk+ΛYjt+δXk=δXk+Ak+1TY+αXkjtxjtkAk+1TYXkTjt+δXk.
With the nonnegativity ofAk+1T,X, and Y,(35)xjtk≥0⇒xjtk+1≥0.
Because(36)xjt←maxxjt−x¯jtATAX+αXXTX+ΛYjt+δ∇XDFXjt,ε,(37)xjtk+1>0.Case 2.
For∇XAk+1,Xkjt<0,(38)xjtk+1=xjtk−maxxijk,δ∇XAk+1,Xkjt∇ADFAij++δXk.
Becausemaxxijk,δ>0 and ∇XAk+1,Xkjt<0, equation (27) holds.
Similarly, the statement ofA can be proven.


### 3.2.1. Convergence Analysis

To analyze the convergence of Algorithm2, we first prove its nonincreasing property under the update rules in equations (25) and (26), i. e.,(39)LAk+1,Xk+1≤LAk,Xk.Because the algorithm solves Problem (1.1) in an alternating method,LAk and LXk can be analyzed separately. Because the update rule of A is similar to that of X, it is sufficient to prove the nonincreasing property of LXk.By making use of the auxiliary function approach [3], the nonincreasing property of LXk can be proved by showing that(40)LXk+1≤GAk+1,Xk+1≤GAk+1,Xk≤GAk,Xk≤LXk,where the auxiliary function G is defined in a similar function by [22](41)GX,Xk=LXk+trX−XkT∇XLXk+12trX−XkTDX−Xk,(42)X=x1x2⋱xT∈ℝ+TJ×J,where xii=1,2,…,T denotes the ith column of X. Meanwhile, we define(43)∇XLXk=∇XLXk1∇XLXk2⋱∇XLXkT∈ℝ+TJ×J,where ∇XLXki denotes the ith column of ∇XLXk, and(44)D=diagD1,D2,…,DT∈ℝ+TJ×TJ,where Dtt=1,…,T denotes a diagonal matrix with its diagonal entries(45)djjt=ATAX+αXXTX+ΛYjt+δXkx¯jtk,ifj∈ln,0,otherwise.Here,ln=j|xjtk>0,∇XLXkjt≠0,orxjtk=0,∇XLXkjt<0 denotes the set of non-KKT indices in the nth column of Xk.Using the Taylor series, the expansion formulation forLX can be expressed as follows:(46)LX=LXk+trX−XkT∇XLXk+12trX−XkTLX−Xk+ξXk,where L=∇X2LXk and ξXk represent the higher components of the Taylor series, and(47)∇X2LXk=∇X2LXk1∇X2LXk2⋱∇X2LXkT∈ℝ+TJ×TJ,where ∇X2LXkn denotes the nth column of ∇X2LXk. Next, for the function GX,Xk, which is an auxiliary of LX, we first have to prove the following Lemma 2. As shown later, Ak and Xk must be bounded for LAk,Xk. The boundedness of Ak and Xk will be proven in Theorem 3.Lemma 2.
Given sufficiently largeδXk and the boundedness of Ak and, then(1)
GX,X=LX.(2)
GXk,Xk=LXk.(3)
GX,X≤GX,Xk. Moreover, GX,X=GX,Xk only if Xk satisfies the KKT conditions.(4)
GX,Xk≤GXk,Xk, and the equality holds only if Xk satisfies the KKT conditions.Proof.
(1) and (2) are obvious from equation (41).
From (3), becauseGX,X=LX, we have(48)GX,Xk−GX,X=GX,X−LX=12trX−XkTD−LX−Xk−ξXk=∑t=1Txt−xtkTDt−Ltx−xtk−ξXk,(49)Dt−Lt=D¯t+δXkD^t−Lt,where D¯t and δXkD^t are diagonal matrices that sum up to Dt with(50)d¯jjt=∇XDFXjt+x¯jtk,ifj∈ln,0,otherwise,d^jjt=1x¯jtk,ifj∈ln,0,otherwise.
Letxt−xtk=x^t; then,(51)GX,Xk−LX=12∑t=1T∑j=1Jx^jt2d¯jjt+δXk∑j=1Jx^jt2d^jjt−12∑t=1Tx^tkTLtx^t−ξXk.
With the boundedness ofAk and Xk and using sufficiently large δXk, the inequality GX,X≤GX,Xk can be guaranteed.
Additionally, it is clear thatGX,X=GX,Xk only if Xk satisfies the KKT conditions demonstrated in equation (48) regardless of ξXk. Because δXk is a variable, GX,X=GX,Xk holds only if X=Xk in equation (51) as a result of the update rule of X and the boundedness of Ak and Xk.
(4) Using(52)GXk,Xk−GX,Xk=−trX−XkT∇XLXkX−Xk−12trX−XkTDX−Xk,and the positive semidefinite of matrix D, GX,Xk is a strictly convex function with respect to Xk; it has a unique minimum that satisfies(53)DX−Xk+∇XLXk=0.
The aforementioned consideration results in the following:(54)GXk,Xk−GX,Xk=12trX−XkTDX−Xk≥0.
IfXk satisfies the KKT conditions, then X=Xk, and therefore, GXk,Xk=GX,Xk.
Conversely, assuming that the equality holds, butXk does not satisfy the KKT conditions, there is at least an index such that(55)xjt≠xjtk,∇XDFXjt++δXkx¯jtk≥δXkx¯jtk.
Furthermore, ifx¯jtk=0, then xjt=xjtk, which contradicts the equality assumption. Therefore, x¯jtk≠0. As a result,(56)GXk,Xk−GX,Xk≥xjt−xjtkδXkx¯jtk>0,which contradicts the assumption.Theorem 1.
Given a sufficiently large numberδXk and supposing that Ak and Xk are bounded, then for the auxiliary function G,(1)
LXk+1≤LXk under the update rule in equation (4) and the equality holds only if Xk satisfies the KKT conditions(2)
LAk+1≤LAk under the update rule in equation (3) and the equality holds only if Ak satisfies the KKT conditions(3)
LAk+1,Xk+1≤LAk,Xk and the equality holds only if Ak,Xk is a stationary pointProof.
(1) and (2) can be obtained directly from Lemma1. (3) By combining (1) and (2), it is obvious that(57)LAk+1,Xk+1≤LAk,Xk,with LAk+1,Xk+1=LAk,Xk only if Ak and Xk satisfy the KKT conditions, i.e., Ak,Xk is a stationary point.Theorem 2.
The objective function (22) is nonincreasing under the multiplicative update rules in equations (25) and (26).Proof.
This theorem is the corollary of Theorem1.Lemma 3.
SequenceAk,Xk,Λk generated by Algorithm 2 is bounded and LAk,Xk is lower bounded.Proof.
Let us verify Assumptions A1–A3 (assumptions A1–A3 can be found in the works of [23]). Assumption A1 holds because of the coercivity of ⋅F. Assumptions A2 and A3 hold because all the coefficient matrices are identity matrices. Additionally, because ⋅F is Lipschitz differentiable, Ak,Xk,Λk generated by Algorithm 3 is bounded, and LAk,Xk is lower bounded according to Theorem 2 presented in the works of [23].
The convergence property of the modified update rules in equations (27) and (29) is stated as follows:Theorem 3.
Given a sufficiently large numberδXk, every limit point of this algorithm is a stationary point of the following optimization problem:(58)minA≥ε,X≥εLA,X=12Y−AXF2+α2XXT−IF2+Λ⋅XXT−I.Proof.
BecauseAk,Xkk=1,…,∞ are in a compact set, the sequence Ak,Xk generated by Algorithm 2 has at least one limit point. Gillis and Glineur [23] introduced a similar analysis that every limit point of this algorithm is a stationary point of the optimization problem (41).


## 3.2.1. Convergence Analysis

To analyze the convergence of Algorithm2, we first prove its nonincreasing property under the update rules in equations (25) and (26), i. e.,(39)LAk+1,Xk+1≤LAk,Xk.Because the algorithm solves Problem (1.1) in an alternating method,LAk and LXk can be analyzed separately. Because the update rule of A is similar to that of X, it is sufficient to prove the nonincreasing property of LXk.By making use of the auxiliary function approach [3], the nonincreasing property of LXk can be proved by showing that(40)LXk+1≤GAk+1,Xk+1≤GAk+1,Xk≤GAk,Xk≤LXk,where the auxiliary function G is defined in a similar function by [22](41)GX,Xk=LXk+trX−XkT∇XLXk+12trX−XkTDX−Xk,(42)X=x1x2⋱xT∈ℝ+TJ×J,where xii=1,2,…,T denotes the ith column of X. Meanwhile, we define(43)∇XLXk=∇XLXk1∇XLXk2⋱∇XLXkT∈ℝ+TJ×J,where ∇XLXki denotes the ith column of ∇XLXk, and(44)D=diagD1,D2,…,DT∈ℝ+TJ×TJ,where Dtt=1,…,T denotes a diagonal matrix with its diagonal entries(45)djjt=ATAX+αXXTX+ΛYjt+δXkx¯jtk,ifj∈ln,0,otherwise.Here,ln=j|xjtk>0,∇XLXkjt≠0,orxjtk=0,∇XLXkjt<0 denotes the set of non-KKT indices in the nth column of Xk.Using the Taylor series, the expansion formulation forLX can be expressed as follows:(46)LX=LXk+trX−XkT∇XLXk+12trX−XkTLX−Xk+ξXk,where L=∇X2LXk and ξXk represent the higher components of the Taylor series, and(47)∇X2LXk=∇X2LXk1∇X2LXk2⋱∇X2LXkT∈ℝ+TJ×TJ,where ∇X2LXkn denotes the nth column of ∇X2LXk. Next, for the function GX,Xk, which is an auxiliary of LX, we first have to prove the following Lemma 2. As shown later, Ak and Xk must be bounded for LAk,Xk. The boundedness of Ak and Xk will be proven in Theorem 3.Lemma 2.
Given sufficiently largeδXk and the boundedness of Ak and, then(1)
GX,X=LX.(2)
GXk,Xk=LXk.(3)
GX,X≤GX,Xk. Moreover, GX,X=GX,Xk only if Xk satisfies the KKT conditions.(4)
GX,Xk≤GXk,Xk, and the equality holds only if Xk satisfies the KKT conditions.Proof.
(1) and (2) are obvious from equation (41).
From (3), becauseGX,X=LX, we have(48)GX,Xk−GX,X=GX,X−LX=12trX−XkTD−LX−Xk−ξXk=∑t=1Txt−xtkTDt−Ltx−xtk−ξXk,(49)Dt−Lt=D¯t+δXkD^t−Lt,where D¯t and δXkD^t are diagonal matrices that sum up to Dt with(50)d¯jjt=∇XDFXjt+x¯jtk,ifj∈ln,0,otherwise,d^jjt=1x¯jtk,ifj∈ln,0,otherwise.
Letxt−xtk=x^t; then,(51)GX,Xk−LX=12∑t=1T∑j=1Jx^jt2d¯jjt+δXk∑j=1Jx^jt2d^jjt−12∑t=1Tx^tkTLtx^t−ξXk.
With the boundedness ofAk and Xk and using sufficiently large δXk, the inequality GX,X≤GX,Xk can be guaranteed.
Additionally, it is clear thatGX,X=GX,Xk only if Xk satisfies the KKT conditions demonstrated in equation (48) regardless of ξXk. Because δXk is a variable, GX,X=GX,Xk holds only if X=Xk in equation (51) as a result of the update rule of X and the boundedness of Ak and Xk.
(4) Using(52)GXk,Xk−GX,Xk=−trX−XkT∇XLXkX−Xk−12trX−XkTDX−Xk,and the positive semidefinite of matrix D, GX,Xk is a strictly convex function with respect to Xk; it has a unique minimum that satisfies(53)DX−Xk+∇XLXk=0.
The aforementioned consideration results in the following:(54)GXk,Xk−GX,Xk=12trX−XkTDX−Xk≥0.
IfXk satisfies the KKT conditions, then X=Xk, and therefore, GXk,Xk=GX,Xk.
Conversely, assuming that the equality holds, butXk does not satisfy the KKT conditions, there is at least an index such that(55)xjt≠xjtk,∇XDFXjt++δXkx¯jtk≥δXkx¯jtk.
Furthermore, ifx¯jtk=0, then xjt=xjtk, which contradicts the equality assumption. Therefore, x¯jtk≠0. As a result,(56)GXk,Xk−GX,Xk≥xjt−xjtkδXkx¯jtk>0,which contradicts the assumption.Theorem 1.
Given a sufficiently large numberδXk and supposing that Ak and Xk are bounded, then for the auxiliary function G,(1)
LXk+1≤LXk under the update rule in equation (4) and the equality holds only if Xk satisfies the KKT conditions(2)
LAk+1≤LAk under the update rule in equation (3) and the equality holds only if Ak satisfies the KKT conditions(3)
LAk+1,Xk+1≤LAk,Xk and the equality holds only if Ak,Xk is a stationary pointProof.
(1) and (2) can be obtained directly from Lemma1. (3) By combining (1) and (2), it is obvious that(57)LAk+1,Xk+1≤LAk,Xk,with LAk+1,Xk+1=LAk,Xk only if Ak and Xk satisfy the KKT conditions, i.e., Ak,Xk is a stationary point.Theorem 2.
The objective function (22) is nonincreasing under the multiplicative update rules in equations (25) and (26).Proof.
This theorem is the corollary of Theorem1.Lemma 3.
SequenceAk,Xk,Λk generated by Algorithm 2 is bounded and LAk,Xk is lower bounded.Proof.
Let us verify Assumptions A1–A3 (assumptions A1–A3 can be found in the works of [23]). Assumption A1 holds because of the coercivity of ⋅F. Assumptions A2 and A3 hold because all the coefficient matrices are identity matrices. Additionally, because ⋅F is Lipschitz differentiable, Ak,Xk,Λk generated by Algorithm 3 is bounded, and LAk,Xk is lower bounded according to Theorem 2 presented in the works of [23].
The convergence property of the modified update rules in equations (27) and (29) is stated as follows:Theorem 3.
Given a sufficiently large numberδXk, every limit point of this algorithm is a stationary point of the following optimization problem:(58)minA≥ε,X≥εLA,X=12Y−AXF2+α2XXT−IF2+Λ⋅XXT−I.Proof.
BecauseAk,Xkk=1,…,∞ are in a compact set, the sequence Ak,Xk generated by Algorithm 2 has at least one limit point. Gillis and Glineur [23] introduced a similar analysis that every limit point of this algorithm is a stationary point of the optimization problem (41).


## 3.3. Auxiliary Variable-Based ADM Extension to ONMF

To facilitate the efficient use of ADM, we first introduce two auxiliary variables,U and V, and consider the following equivalent model:(59)minX,Y,U,V12Y−AXF2s.t.A=U,X=V,U≥0,V≥0,XXT=I,where A,U∈ℝm×q and Y,V∈ℝq×n. The augmented Lagrangian of equation (59) is as follows:(60)LAA,X,U,V,Λ,Π=12Y−AXF2+α2XXT−IF2+Λ⋅A−U+Π⋅X−V+ρ12A−UF2+ρ22X−VF2,where Λ∈ℝm×q,Π∈ℝq×n are Lagrangian multipliers, α,ρ1,ρ2>0 are penalty parameters, and A⋅B≔∑i,jaijbij for matrices A and B of the same size.The alternating direction method for equation (59) is derived by successively minimizing LA with respect to A,X,U,V, one at a time, while fixing others at their more recent values, i. e.,(61)Ak+1=argminLAA,Xk,Uk,Vk,Λk,Πk,Xk+1=argminLAAk+1,X,Uk,Vk,Λk,Πk,Uk+1=argminU≥0LAAk+1,Xk+1,U,Vk,Λk,Πk,Vk+1=argminU≥0LAAk+1,Xk+1,Uk+1,V,Λk,Πk,and then updating the multipliers Λ and Π. The introduction of the two auxiliary variables, U and V, makes it easy to carry out each of the alternating minimization steps. Specifically, these steps can be expressed in closed form as follows:(62)aij←aijYXT+ρ1UijAXXT+Π+ρ1Aij+δω,(63)xjt←xjtATY+αX+ρ2VjtATAX+αXXTX+Λ+ρ2Xjt+δω,(64)U←P+A+Λρ1,(65)V←P+X+Πρ2,(66)Λ←Λ+γρ1A−U,(67)Π←Π+γρ2X−V.The detailed procedure for the algorithm of Problem59 can be seen in Algorithm 3.


## 4. Experiments and Results

In this section, we analyze and compare Algorithms2 and 3 with Algorithms 1 and 4 (ONP-MF [18]). All experiments are performed using the Windows 7 operating system and MATLAB v8.1 (R2013a) running on a Lenovo laptop with an Intel Core(TM)i5-3470 CPU at 3.2 GHz and 4 GB of memory.


### 4.1. Test on Images

To visualize the results, we use examples of ONMF arising in image processing to analyze the quality of the solutions. Two grayscale images are tested using Algorithms1–3 and 4. The two grayscale images, a panda in Figure 1(a) and a cat in Figure 1(b), have resolutions of 1200×1920 and 128×128, respectively.Figure 1
Original images. (a) Panda. (b) Cat.
(a)(b)We evaluate the degree of approximation and orthogonality using two indices:(68)Relative residual errorRRE:Y−AXFYF,Orthogonality residual errorORE:XXT−IF.We settol=10−5, δ=10−8, and α=0.1, and maxiter = 500 for the image of the cat and maxiter = 100 for the image of the panda, after which we repeat the process 10 times and calculate the relative errors and orthogonality. The results are given in Tables 1–5, and the recovered images are shown in Figures 2 and 3.Table 1
Comparison ofREE testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 4100.2256 ± 0.02160.2096±0.02760.2413 ± 0.03550.2191 ± 0.0227200.1973 ± 0.01730.1677 ± 0.01830.1811 ± 0.01790.1761 ± 0.0166400.1876 ± 0.01250.1310 ± 0.01560.1665 ± 0.01420.1296±0.0114600.1722 ± 0.00570.1278 ± 0.00470.1542 ± 0.01130.1204±0.0086800.1580 ± 0.00630.1105 ± 0.00810.1401 ± 0.00860.1079±0.0068Table 2
Comparison ofREE testing on image cat.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 450.1906 ± 0.03140.1855±0.01910.1957 ± 0.02410.1947 ± 0.0264100.1761 ± 0.02450.1638±0.01760.1740 ± 0.02110.1711 ± 0.0216150.1659 ± 0.01740.1513 ± 0.01520.1601 ± 0.01870.1473 ± 0.0182200.1584 ± 0.01120.1362±0.01210.1474 ± 0.01580.1384 ± 0.0155250.1521 ± 0.01210.1252 ± 0.01020.1392 ± 0.01290.1234±0.0133300.1477 ± 0.01890.1175 ± 0.00930.1292 ± 0.01110.1103±0.0115Table 3
Comparison ofORE testing on image cat.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 451.1121 ± 0.32170.4588±0.22110.9992 ± 0.00010.9091 ± 0.3015102.7271 ± 0.24730.5600±0.20080.9994 ± 0.00020.9524 ± 0.2182154.3523 ± 0.18200.5681±0.16560.9996 ± 0.00010.9677 ± 0.1796205.8419 ± 0.15560.5727±0.13570.9996 ± 0.00030.9756 ± 0.1562257.8002 ± 0.14820.5799±0.12290.9997 ± 0.00020.9804 ± 0.1400309.4566 ± 0.13310.5860±0.09800.9998 ± 0.00020.9836 ± 0.1280Table 4
Comparison ofORE testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 410372.6551 ± 10.25460.5457±0.17550.9955 ± 0.12570.9210 ± 0.154620686.1020 ± 40.11200.5601±0.18150.9979 ± 0.15220.9664 ± 0.1753401.3621×103 ± 70.18850.6573±0.18000.9995 ± 0.25610.9718 ± 0.2051601.9141×103 ± 114.28500.6607±0.21110.9999 ± 0.21010.9957 ± 0.2689802.6000×103 ± 160.47820.8610±0.15661.0000 ± 0.20080.9980 ± 0.1280Table 5
Comparison of cpu time (seconds) testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 41093.1131 ± 5.2421427.2460 ± 10.854187.0101 ± 6.010515.7797±2.30852096.2495 ± 4.2680470.1573 ± 15.963098.3510 ± 8.150720.9627±3.03434099.1576 ± 8.0485526.0207 ± 18.0012105.1883 ± 9.205736.4698±5.224060101.6872 ± 10.5701579.2561 ± 23.2007112.1097 ± 10.000558.2959±8.417880104.0500 ± 9.4910819.0018 ± 30.0804120.2542 ± 13.507690.9566±13.0304Figure 2
Recovered128×128 cat by Algorithms 1–3, and Algorithm 4. (a) J = 5. (b) J = 5. (c) J = 5. (d) J = 5. (e) J = 15. (f) J = 15. (g) J = 15. (h) J = 15. (i) J = 30. (j) J = 30. (k) J = 30. (l) J = 30.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 3
Recovered1200×1920 panda by Algorithms 1–3, and Algorithm 4. (a) J = 10. (b) J = 10. (c) J = 10. (d) J = 10. (e) J = 20. (f) J = 20. (g) J = 20. (h) J = 20. (i) J = 40. (j) J = 40. (k) J = 40. (l) J = 40. (m) J = 80. (n) J = 80. (o) J = 80. (p) J = 80.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Tables1 and 2 show that Algorithm 2 performs slightly better than Algorithms 1 and 3 in terms of relative errors in most cases. Moreover, we observe that Algorithm 4 converges in least iterations. From Tables 3 and 4, we can see that Algorithm 2 demonstrates the best performance compared with Algorithms 1, 3 and 4 in terms of orthogonality. From Table 5, we observe that Algorithm 4 requires the lowest computational cost compared with the rest of algorithms. From Figures 2 and 3, we observe that Algorithm 2 achieves a better recovery quality of images compared with Algorithms 1 and 3, and the recovery quality of Algorithm 4 is the best.


### 4.2. Clustering Capability

We test the clustering capability of the three ONMF algorithms on the ALLAML gene-sample data provided by [4] and compare the performance of the algorithms. To measure the clustering performance, we use purity and entropy as our performance metrics. We expect that these metrics will provide us with good insights on how our algorithm works.Entropy measures how classes of genes are distributed within each cluster [25]. It is defined as follows:(69)entropy≡1Nlog2S∑r=1R∑s=1Scrslog2crscr,where N denotes the number of samples, crs denotes the number of samples in the rth cluster that belong to the sth class, and cr denotes the size of the rth cluster.Purity is the most commonly used metric. It measures the extent to which each cluster contains data points from primarily one class [25]. Therefore, the higher the value of purity, the better the clustering performance. Purity is defined as follows:(70)purity≡1N∑r=1Rmaxscrs.Algorithms2 and 3 are compared with Algorithms 1 and 4, and the comparisons are shown in Tables 6 and 7. From Tables 6 and 7, we observe that our Algorithms 2 and 3 clustering achieves better purity and entropy results on the leukemia data sets compared with the performance of Algorithms 1 and 4. There exist slight differences in the relative performance of purity and entropy in our comparison, that is, higher purity values do not necessarily correspond to lower entropy values. This is because entropy takes into account the entire distribution of the genes in a particular cluster and not only the largest class as in the computation of purity. In summary, the comparison shows that Algorithms 2 and 3 are viable and competitive in clustering tasks.Table 6
Performancek = 2 on the leukemia data matrix.
Algorithm1Algorithm2Algorithm3Algorithm 4Purity0.9216± 0.12560.9501±0.07610.9266± 0.11270.6708± 0.0950Entropy0.3053± 0.14730.2390±0.06550.3317± 0.09520.6653± 0.0958Table 7
Performancek = 3 on the leukemia data matrix.
Algorithm1Algorithm2Algorithm3Algorithm 4Purity0.8912± 0.10240.9360±0.15860.9218± 0.18950.8256± 0.1179Entropy0.2500± 0.05670.1796±0.04650.2054± 0.07230.4132± 0.0590


## 4.1. Test on Images

To visualize the results, we use examples of ONMF arising in image processing to analyze the quality of the solutions. Two grayscale images are tested using Algorithms1–3 and 4. The two grayscale images, a panda in Figure 1(a) and a cat in Figure 1(b), have resolutions of 1200×1920 and 128×128, respectively.Figure 1
Original images. (a) Panda. (b) Cat.
(a)(b)We evaluate the degree of approximation and orthogonality using two indices:(68)Relative residual errorRRE:Y−AXFYF,Orthogonality residual errorORE:XXT−IF.We settol=10−5, δ=10−8, and α=0.1, and maxiter = 500 for the image of the cat and maxiter = 100 for the image of the panda, after which we repeat the process 10 times and calculate the relative errors and orthogonality. The results are given in Tables 1–5, and the recovered images are shown in Figures 2 and 3.Table 1
Comparison ofREE testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 4100.2256 ± 0.02160.2096±0.02760.2413 ± 0.03550.2191 ± 0.0227200.1973 ± 0.01730.1677 ± 0.01830.1811 ± 0.01790.1761 ± 0.0166400.1876 ± 0.01250.1310 ± 0.01560.1665 ± 0.01420.1296±0.0114600.1722 ± 0.00570.1278 ± 0.00470.1542 ± 0.01130.1204±0.0086800.1580 ± 0.00630.1105 ± 0.00810.1401 ± 0.00860.1079±0.0068Table 2
Comparison ofREE testing on image cat.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 450.1906 ± 0.03140.1855±0.01910.1957 ± 0.02410.1947 ± 0.0264100.1761 ± 0.02450.1638±0.01760.1740 ± 0.02110.1711 ± 0.0216150.1659 ± 0.01740.1513 ± 0.01520.1601 ± 0.01870.1473 ± 0.0182200.1584 ± 0.01120.1362±0.01210.1474 ± 0.01580.1384 ± 0.0155250.1521 ± 0.01210.1252 ± 0.01020.1392 ± 0.01290.1234±0.0133300.1477 ± 0.01890.1175 ± 0.00930.1292 ± 0.01110.1103±0.0115Table 3
Comparison ofORE testing on image cat.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 451.1121 ± 0.32170.4588±0.22110.9992 ± 0.00010.9091 ± 0.3015102.7271 ± 0.24730.5600±0.20080.9994 ± 0.00020.9524 ± 0.2182154.3523 ± 0.18200.5681±0.16560.9996 ± 0.00010.9677 ± 0.1796205.8419 ± 0.15560.5727±0.13570.9996 ± 0.00030.9756 ± 0.1562257.8002 ± 0.14820.5799±0.12290.9997 ± 0.00020.9804 ± 0.1400309.4566 ± 0.13310.5860±0.09800.9998 ± 0.00020.9836 ± 0.1280Table 4
Comparison ofORE testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 410372.6551 ± 10.25460.5457±0.17550.9955 ± 0.12570.9210 ± 0.154620686.1020 ± 40.11200.5601±0.18150.9979 ± 0.15220.9664 ± 0.1753401.3621×103 ± 70.18850.6573±0.18000.9995 ± 0.25610.9718 ± 0.2051601.9141×103 ± 114.28500.6607±0.21110.9999 ± 0.21010.9957 ± 0.2689802.6000×103 ± 160.47820.8610±0.15661.0000 ± 0.20080.9980 ± 0.1280Table 5
Comparison of cpu time (seconds) testing on image panda.
RankJAlgorithm1Algorithm2Algorithm3Algorithm 41093.1131 ± 5.2421427.2460 ± 10.854187.0101 ± 6.010515.7797±2.30852096.2495 ± 4.2680470.1573 ± 15.963098.3510 ± 8.150720.9627±3.03434099.1576 ± 8.0485526.0207 ± 18.0012105.1883 ± 9.205736.4698±5.224060101.6872 ± 10.5701579.2561 ± 23.2007112.1097 ± 10.000558.2959±8.417880104.0500 ± 9.4910819.0018 ± 30.0804120.2542 ± 13.507690.9566±13.0304Figure 2
Recovered128×128 cat by Algorithms 1–3, and Algorithm 4. (a) J = 5. (b) J = 5. (c) J = 5. (d) J = 5. (e) J = 15. (f) J = 15. (g) J = 15. (h) J = 15. (i) J = 30. (j) J = 30. (k) J = 30. (l) J = 30.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 3
Recovered1200×1920 panda by Algorithms 1–3, and Algorithm 4. (a) J = 10. (b) J = 10. (c) J = 10. (d) J = 10. (e) J = 20. (f) J = 20. (g) J = 20. (h) J = 20. (i) J = 40. (j) J = 40. (k) J = 40. (l) J = 40. (m) J = 80. (n) J = 80. (o) J = 80. (p) J = 80.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Tables1 and 2 show that Algorithm 2 performs slightly better than Algorithms 1 and 3 in terms of relative errors in most cases. Moreover, we observe that Algorithm 4 converges in least iterations. From Tables 3 and 4, we can see that Algorithm 2 demonstrates the best performance compared with Algorithms 1, 3 and 4 in terms of orthogonality. From Table 5, we observe that Algorithm 4 requires the lowest computational cost compared with the rest of algorithms. From Figures 2 and 3, we observe that Algorithm 2 achieves a better recovery quality of images compared with Algorithms 1 and 3, and the recovery quality of Algorithm 4 is the best.


## 4.2. Clustering Capability

We test the clustering capability of the three ONMF algorithms on the ALLAML gene-sample data provided by [4] and compare the performance of the algorithms. To measure the clustering performance, we use purity and entropy as our performance metrics. We expect that these metrics will provide us with good insights on how our algorithm works.Entropy measures how classes of genes are distributed within each cluster [25]. It is defined as follows:(69)entropy≡1Nlog2S∑r=1R∑s=1Scrslog2crscr,where N denotes the number of samples, crs denotes the number of samples in the rth cluster that belong to the sth class, and cr denotes the size of the rth cluster.Purity is the most commonly used metric. It measures the extent to which each cluster contains data points from primarily one class [25]. Therefore, the higher the value of purity, the better the clustering performance. Purity is defined as follows:(70)purity≡1N∑r=1Rmaxscrs.Algorithms2 and 3 are compared with Algorithms 1 and 4, and the comparisons are shown in Tables 6 and 7. From Tables 6 and 7, we observe that our Algorithms 2 and 3 clustering achieves better purity and entropy results on the leukemia data sets compared with the performance of Algorithms 1 and 4. There exist slight differences in the relative performance of purity and entropy in our comparison, that is, higher purity values do not necessarily correspond to lower entropy values. This is because entropy takes into account the entire distribution of the genes in a particular cluster and not only the largest class as in the computation of purity. In summary, the comparison shows that Algorithms 2 and 3 are viable and competitive in clustering tasks.Table 6
Performancek = 2 on the leukemia data matrix.
Algorithm1Algorithm2Algorithm3Algorithm 4Purity0.9216± 0.12560.9501±0.07610.9266± 0.11270.6708± 0.0950Entropy0.3053± 0.14730.2390±0.06550.3317± 0.09520.6653± 0.0958Table 7
Performancek = 3 on the leukemia data matrix.
Algorithm1Algorithm2Algorithm3Algorithm 4Purity0.8912± 0.10240.9360±0.15860.9218± 0.18950.8256± 0.1179Entropy0.2500± 0.05670.1796±0.04650.2054± 0.07230.4132± 0.0590


## 5. Conclusion

In this study, we develop two efficient modified algorithms based on the ADM algorithms for the ONMF problem. We then prove that any sequence of solutions generated by the modified updates of Algorithms2 and 3 has at least one convergent subsequence, and the limit of any convergent subsequence is a stationary point of the corresponding optimization problem. From the numerical results, a notable observation is that our Algorithms 2 and 3 are feasible and efficient, especially in the properties of the solution, that is, the orthogonality.The framework for the convergence analysis of the updates for Algorithm3 presented in this article may be applied to the ADM algorithms for NMF [26].



---

*Source: 8490147-2021-07-01.xml*