{'fraction_non_alphanumeric': 0.046485088253195375, 'fraction_numerical': 0.011817812943801988, 'mean_word_length': 4.210491543340381, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonnegative Factorization and The Maximum Edge Biclique Problem

Nonnegative Matrix Factorization (NMF) is a data analysis technique which allows compression and interpretation of nonnegative data. NMF became widely studied after the publication of the seminal paper by Lee and Seung (Learning the Parts of Objects by Nonnegative Matrix Factorization, Nature, 1999, vol. 401, pp. 788--791), which introduced an algorithm based on Multiplicative Updates (MU). More recently, another class of methods called Hierarchical Alternating Least Squares (HALS) was introduced that seems to be much more efficient in practice. In this paper, we consider the problem of approximating a not necessarily nonnegative matrix with the product of two nonnegative matrices, which we refer to as Nonnegative Factorization (NF); this is the subproblem that HALS methods implicitly try to solve at each iteration. We prove that NF is NP-hard for any fixed factorization rank, using a reduction to the maximum edge biclique problem. We also generalize the multiplicative updates to NF, which allows us to shed some light on the differences between the MU and HALS algorithms for NMF and give an explanation for the better performance of HALS. Finally, we link stationary points of NF with feasible solutions of the biclique problem to obtain a new type of biclique finding algorithm (based on MU) whose iterations have an algorithmic complexity proportional to the number of edges in the graph, and show that it performs better than comparable existing methods.


Introduction
(Approximate) Nonnegative Matrix Factorization (NMF) is the problem of approximating a given nonnegative matrix by the product of two low-rank nonnegative matrices: given a matrix M ≥ 0, one has to compute two low-rank matrices V, W ≥ 0 such that M ≈ V W . (1.1) This problem was first introduced in 1994 by Paatero and Tapper [25], and more recently received a considerable interest after the publication of two papers by Lee and Seung [21,22]. It is now well established that NMF is useful in the framework of compression and interpretation of nonnegative data ; it has for example been applied in analysis of image databases, text mining, interpretation of spectra, computational biology and many other applications (see e.g. [2,9,11] and references therein).
How can one interpret the outcome of a NMF? Assume each column M :j of matrix M represents an element of a data set: expression (1.1) can be equivalently written as where each element M :j is decomposed into a nonnegative linear combination (with weights W kj ) of nonnegative basis elements ({V :k }, the columns of V ). Nonnegativity of V allows interpretation of the basis elements in the same way as the original nonnegative elements in M , which is crucial in applications where the nonnegativity property is a requirement (e.g. where elements are images described by pixel intensities or texts represented by vectors of word counts). Moreover, nonnegativity of the weight matrix W corresponds to an essentially additive reconstruction which leads to a partbased representation: basis elements will represent similar parts of the columns of M . Sparsity is another important consideration: finding sparse factors improves compression and leads to a better part-based representation of the data [19]. We start this paper with a brief introduction to the NMF problem: Section 2 recalls existing complexity results, introduces two well-known classes of methods: multiplicative updates [22] and hierarchical alternating least squares [6] and proposes a simple modification to guarantee their convergence. The central problem studied in this paper, Nonnegative Factorization (NF), is a generalization of NMF where the matrix to be approximated with the product of two low-rank nonnegative matrices is not necessarily nonnegative. NF is introduced in Section 3, where it is shown to be NP-hard for any given factorization rank, using a reduction to the problem of finding a maximum edge biclique. Stationary points of the NF problem used in that reduction are also studied. This section ends with a generalization of the NMF multiplicative updates rules to the NF problem and a proof of their convergence. This allows us to shed new light on the standard multiplicative updates for NMF: a new interpretation is given in Section 4, which explains the relatively poor performance of these methods and hints at possible improvements. Finally, Section 5 introduces a new type of biclique finding algorithm that relies on the application of multiplicative updates to the equivalent NF problem considered earlier. This algorithm only requires a number of operations proportional to the number of edges of the graph per iteration, and is shown to perform well when compared to existing methods.

Nonnegative Matrix Factorization (NMF)
Given a matrix M ∈ R m×n + and an integer r ∈ N 0 , the NMF optimization problem using the Frobenius norm is defined as R m×n denotes the set of real matrices of dimension m × n; R m×n + the set of nonnegative matrices i.e. R m×n with every entry nonnegative, and 0 the zero matrix of appropriate dimensions. A wide range of algorithms have been proposed to find approximate solutions for this problem (see e.g. [2,4,6,7,10,11,13,24]). Most of them use the fact that although problem (NMF) is not convex, its objective function is convex separately in each of the two factors V and W (which implies that finding the optimal factor V corresponding to a fixed factor W reduces to a convex optimization problem, and vice-versa), and try to find good approximate solutions by using alternating minimization schemes. For instance, Nonnegative Least Squares (NNLS) algorithms can be used to minimize (exactly) the cost function alternatively over factors V and W (see e.g. [5,20]).
Actually, there exist other partitions of the variables that preserve convexity of the alternating minimization subproblems: since the cost function can be rewritten as ||M − r i=1 V :i W i: || F , it is clearly convex as long as variables do not include simultaneously an element of a column of V and an element of the corresponding row of W (i.e. V ki and W il for the same index i). Therefore, given a subset of indexes K ⊆ R = {1, 2, . . . , r}, (NMF) is clearly convex for both the following subsets of variables However, the convexity is lost as soon as one column of V (V :i ) and the corresponding row of W (W i: ) are optimized simultaneously, so that the corresponding minimization subproblem can no longer be efficiently solved up to global optimality.

Complexity
Vavasis studies in [30] the algorithmic complexity of the NMF optimization problem; more specifically, he proves that the following problem, called Exact Nonnegative Matrix Factorization 1 , is NP-hard: (Exact NMF) Given a nonnegative matrix M ≥ 0 of rank k, find, if possible, two nonnegative factors V ≥ 0 and W ≥ 0 of rank k such that M = V W .
The NMF optimization problem is therefore also NP-hard, since when the rank r is equal to the rank k of the matrix M , any optimal solution to the NMF optimization problem can be used to answer the Exact NMF problem (the answer being positive if and only if the optimal objective value of the NMF optimization problem is equal to zero).
The NP-hardness proof for exact NMF relies on its equivalence with a NP-hard problem in polyhedral combinatorics, and requires both the dimensions of matrix M and its rank k to increase to obtain NP-hardness. In contrast, in the special cases when rank k is equal to 1 or 2, the exact NMF problem can always be answered in the affirmative: 1. When k = 1, it is obvious that for any nonnegative rank-one matrix M ≥ 0 there is nonnegative factors v ≥ 0 and w ≥ 0 such that M = vw T .
Moreover, the NMF optimization problem with r = 1 can be solved in polynomial time: the Perron-Frobenius theorem implies that the dominant left and right singular vectors of a nonnegative matrix M are nonnegative, while the Eckart-Young theorem states that the outer product of these dominant singular vectors is the best rank-one approximation of M ; these vectors can be computed in polynomial-time using for example the singular value decomposition [15].
1 This is closely related to the nonnegative rank of matrix M , which is the minimum value of r for which there exists V ∈ R m×r + and W ∈ R r×n + such that M = V W (see [1]).
2. When nonnegative matrix M has rank 2, Thomas has shown [29] that exact NMF is also always possible (see also [8]). The fact that any rank-two nonnegative matrix can be exactly factorized as the product of two rank-two nonnegative matrices can be explained geometrically as follows: viewing columns of M as points in R m , the fact that M has rank 2 implies that the set of its columns belongs to a two-dimensional subspace. Furthermore, because these columns are nonnegative, they belong to a two-dimensional pointed cone, see Figure 1. Since such a cone is always spanned by two extremes vectors, this implies that all columns of M can be represented exactly as nonnegative linear combinations of two nonnegative vectors, and therefore the exact NMF is always possible 2 . Moreover, these two extreme columns can easily be computed in polynomial time (using for example the fact that they define an angle of maximum amplitude among all pairs of columns). Hence, when the optimal rank-two approximation of matrix M is nonnegative, the NMF optimization problem with r = 2 can be solved in polynomial time. However, this optimal rank-two approximation is not always nonnegative, so that the complexity of the NMF optimization in the case r = 2 is not known. Furthermore, to the best of our knowledge, the complexity of the exact NMF problem and the NMF optimization problem are still unknown for any fixed rank r or k greater than 3.

Multiplicative Updates (MU)
In their seminal paper [22], Lee and Seung propose multiplicative update rules that aim at minimizing the Frobenius norm between M and V W . To understand the origin of these rules, consider the Karush-Kuhn-Tucker first-order optimality conditions for (NMF) where • is the Hadamard (component-wise) product between two matrices, and From these equalities, Lee and Seung derive the following simple multiplicative update rules (where [.] [.] is Hadamard (component-wise) division) for which they are able to prove a monotonicity property: ). The Frobenius norm ||M −V W || F is nonincreasing under the multiplicative update rules (2.7).
The algorithm based on the alternated application of these rules is not guaranteed to converge to a first-order stationary point (see e.g. [2], and references therein), but a slight modification proposed in [23] achieves this property (roughly speaking, MU is recast as a variable metric steepest descent method and the step length is modified accordingly). We propose another possibility to overcome this problem by replacing the above updates by the following: for any (V, W ) ≥ . Moreover, every limit point of this algorithm is a stationary point of the following optimization problem min Proof. See Section 3.3 of this paper, where the more general Theorem 9 is proved.

Hierarchical Alternating Least Squares (HALS)
Cichoki et al. [6] and independently several other authors [14,18] have proposed to solve the problem of Nonnegative Matrix Factorization by considering successively each rank-one factor V :k W k: while keeping the rest of the variables fixed, which can be expressed as where matrix R k is called the k th residual matrix. Ideally, one would like to find an optimal rank-one factor V :k W k: according to the Frobenius norm, i.e. solve the following problem but, instead of solving this problem directly, these authors propose to optimize column V :k and row W k: separately in an alternating scheme, because the optimal solution to these two (convex) subproblems can be easily computed in closed form, see e.g. [16]: ||W k: || 2 2 (2.12) (2.13) This scheme, which amounts to a block coordinate descent method (for which any cyclic order on the columns of V and the rows of W is admissible), is called Hierarchical Alternating Least Squares (HALS) 3 and it has been observed to work remarkably well in practice: it outperforms, in most cases, the other algorithms for NMF [7,14,16]. Indeed, it combines a low computational cost per iteration (the same as the multiplicative updates) with a relatively fast convergence (significantly faster than the multiplicative updates), see Figure 3 for an example. We will explain later in Section 4 why this algorithm performs much better than the one of Lee and Seung. A potential issue with this method is that, in the course of the optimization process, one of the vectors V :k (or W k: ) and the corresponding rank-one factor V :k W k: may become equal to zero (this happens for example if one of the residuals R k is nonpositive). This then leads to numerical instabilities (the next update is not well-defined) and a rank-deficient approximation (with a rank lower than r). A possible way to overcome this problem is to replace the zero lower bounds on V :k and W k: in (2.12) and (2.13) by a small positive constant, say ≪ 1 (as for the MU), and consider the following subproblems which lead to the modified closed-form update rules: This idea was already suggested in [6] in order to avoid numerical instabilities. In fact, this variant of the algorithm is now well-defined in all situations because (2.14) guarantees V :k > 0 and W k: > 0 at each iteration. Furthermore, one can now easily prove that it converges to a stationary point.
Theorem 3. For every constant > 0, the limit points of the block coordinate descent algorithm initialized with positive matrices and applied to the optimization problem (2.9) are stationary points.
Proof. We use the following result of Powell [27] (see also [3, p.268]): the limit points of the iterates of a block coordinate descent algorithm are stationary points provided that the following two conditions hold: • each block of variables is required to belong to a closed convex set, • the minimum computed at each iteration for a given block of variables is uniquely attained.
The first condition is clearly satisfied here, since V :k and W k: belong respectively to ([ , +∞[) m and ([ , +∞[) n , which are closed convex sets. The second condition holds because subproblems (2.14) can be shown to be strictly convex, so that their optimal value is uniquely attained by the solutions provided by rules (2.15). Strict convexity is due to the fact that the objective function of these problems are sums of quadratic terms, each involving a single variable and having a strictly positive coefficient.

Nonnegative Factorization (NF)
Looking back at subproblem (2.11), i.e. approximating the residual R k with a rank-one term V :k W k: , we have seen that the optimal solution separately for both V :k and W k: can be written in a closed form.
In the previous section, subproblem (2.11) was then solved by a block coordinate descent algorithm. A question arises: Is it possible to do better? i.e. Is it possible to efficiently solve the problem for both vectors simultaneously? In order to answer this question, we introduce the problem of Nonnegative Factorization 4 which is exactly the same as Nonnegative Matrix Factorization except that the matrix to factorize can be any real matrix, i.e. is not necessarily nonnegative. Given M ∈ R m×n and r ∈ N 0 , the Nonnegative Factorization optimization problem using the Frobenius norm is: Of course, this problem is a generalization of (NMF) and is NP-hard as well. However Nonnegative Factorization will be shown below to be NP-hard for any fixed factorization rank (even r = 1), which is not the case of (NMF) (cf. Section 2.1). The proof is based on the reduction to the maximum edge biclique problem.

Complexity
The main result of this section is the NP-hardness result for Nonnegative Factorization for any fixed factorization rank. We first show how the optimization version of the maximum edge biclique problem (MBP) can be formulated as a rank-one Nonnegative Factorization problem (NF-1d). Since the decision version of (MBP) is NP-complete [26], this implies that (NF-1d) is NP-hard. We then prove that (NF) is NP-hard as well using a simple construction.

The Maximum Edge Biclique Problem in Bipartite Graphs
A bipartite graph G b is a graph whose vertices can be divided into two disjoint sets V 1 and V 2 such that there is no edge between two vertices in the same set A biclique K b is a complete bipartite graph i.e. a bipartite graph where all the vertices are connected Finally, the so-called maximum edge biclique problem in a bipartite graph maximizing the number of edges. The decision problem: Given B, does G b contain a biclique with at least B edges? has been shown to be NP-complete [26]. Therefore the corresponding optimization problem is at least NP-hard.
In order to avoid trivialities, we will suppose that each vertex of the graph is connected to at least one other vertex i.e. M b (i, : We denote by |E| the cardinality of E i.e. the number of edges in G b ; note that |E| = ||M b || 2 F . The set of zero values will be called Z = {(i, j) | M b (i, j) = 0}, and its cardinality |Z| satisfies |E| + |Z| = mn.
With this notation, the maximum biclique problem in G b can be formulated as In fact, one can check easily that this objective is equivalent to max v,w ij v i w j since M b , v and w are binary: instead of maximizing the number of edges inside the biclique, one minimizes the number of edges outside. Feasible solutions of (MBP) correspond to bicliques of G b . We will be particularly interested in maximal bicliques. A maximal biclique is a biclique which is not contained in a larger biclique: it is a locally optimal solution of (MBP).
The corresponding rank-one Nonnegative Factorization problem is defined as with the matrix M d defined as To prove NP-hardness of (NF-1d), we are going to show that, if d is sufficiently large, optimal solutions of (NF-1d) coincide with optimal solutions of the corresponding biclique problem (MBP).
By the Eckart-Young theorem, where σ max (M ) is the maximum singular value of M . Clearly, We restate here a well-known result concerning low-rank approximations (see e.g. [16, p. 29]). Lemma 2. The local minima of the best rank-one approximation problem with respect to the Frobenius norm are global minima.
We can now state the main result about the equivalence of (NF-1d) and (MBP).
Theorem 4. For d ≥ |E|, any optimal solution (v,w) of (NF-1d) coincides with an optimal solution of (MBP), i.e. vw is binary and vw ≤ M b .
Proof. We focus on the entries of vw which are positive and define is positive and is an optimal solution of (NF-1d) for M d , (v , w ) is a local minimum of the unconstrained problem i.e. the problem of best rankone approximation. By Lemma 2, this must be a global minimum. This is a contradiction with Lemma 1: (v , w ) should contain at least one nonpositive entry. Therefore, M d = 1 |K|×|L| which implies v w = M d by optimality and then vw is binary and vw ≤ M b .
Intuitively, the reason (NF-1d) is NP-hard is that if one of the −d entries of M d is approximated by a positive value, say p, the corresponding error is d 2 + 2pd + p 2 . Therefore, the larger d, the more expensive it is to approximate −d by a positive number. Because of that, when d increases, negatives values of M d will be approximated by smaller values and eventually by zeros.
Therefore, for each negative entry of M , one has to decide whether to approximate it with zero or with a positive value. Moreover, when a value is approximated by zero, one has to choose which entries of V and W will be equal to zero, as in the biclique problem. We suspect that the hardness of Nonnegative Factorization lies in these combinatorial choices.
We can now answer our initial question: Would it be possible to solve efficiently the problem simultaneously for both vectors (V :k , W k: )? Obviously, unless P=NP, we won't be able to find a polynomial-time algorithm to solve this problem. Therefore, it seems hopeless to improve the HALS algorithm using this approach. Remark 1. Corollary 1 suggests that NMF is a difficult problem for any fixed r ≥ 2. Indeed, even if one was given the optimal solution of an NMF problem except for one rank-one factor, it is not guaranteed that one would be able to find this last factor in polynomial-time since the corresponding residue is not necessarily nonnegative.
We now generalize Theorem 1 to factorizations of arbitrary rank. Proof. Let M b ∈ {0, 1} m×n be the adjacency matrix of a bipartite graph G b and r ≥ 1 the factorization rank of (NF). We define the matrix A b as which is the adjacency matrix of another bipartite graph G b which is nothing but the graph G b repeated r times. A d is defined in the same way as M d , i.e.
Let (V, W ) be the optimal rank-r nonnegative factorization of A d and consider each rank-one factor V :k W k: ≈ R k = A d − i =k V :i W i: : each of them must clearly be an optimal rank-one nonnegative factorization of R k . Since R k ≤ A d , and Lemma 1 holds. Using exactly the same arguments as in Theorem 4, one can show that, ∀k, Therefore, the positive entries of each rank-one factor will correspond to a biclique of G b . By optimality of (V, W ), each rank-one factor must correspond to a maximum biclique of G b since G b is the graph G b repeated r times. Thus (NF-1d) is NP-hard implies (NF) is NP-hard.

Stationary points of (NF-1d)
We have shown that optimal solutions of (NF-1d) coincide with optimal solutions of (MBP) for d ≥ √ E, which are NP-hard to find. In this section, we focus on stationary points of (NF-1d) instead: we show how they are related to the feasible solutions of (MBP). This result will be used in Section 5 to design a new type of biclique finding algorithm.

Definitions and Notation
The KKT conditions of (NF-1d), which define the stationary points, are exactly the same as for (NMF): (v, w) is a stationary point of (NF-1d) if and only if Of course, we are especially interested in nontrivial solutions and we then assume v, w = 0 so that one can check that Given d, we define three sets of rank-one matrices: S d , corresponding to the nontrivial stationary points of (NF-1d), with F , corresponding to the feasible solutions of (MBP), with is a feasible of (MBP)}, and B, corresponding to the maximal bicliques of (MBP), i.e. vw ∈ B if and only if vw ∈ F and vw coincides with a maximal biclique.

Stationarity of Maximal Bicliques
The next theorem states that, for d sufficiently large, the only nontrivial feasible solutions of (MBP) that are stationary points of (NF-1d) are the maximal bicliques.
Proof. vw ∈ B if and only if vw ∈ F and is maximal i.e.
(1) i such that v i = 0 and M d (i, j) = 1, ∀j s.t. w j = 0, Since vw is binary and v, w = 0, the nonzero entries of v and w must be equal to each other. Moreover, This is equivalent to the stationarity conditions for v = 0, cf. (3.6). By symmetry, (2) is equivalent to the stationarity conditions for w.
Theorem 6 implies that, for d sufficiently large, B ⊂ S d . It would be interesting to have the opposite affirmation: for d sufficiently large, any stationary point of (NF-1d) corresponds to a maximal biclique of (MBP). Unfortunately, we will see later that this property does not hold.

Limit points of S d
However, as d goes to infinity, we are going to show that the points in S d get closer to feasible solutions of (MBP). Lemma 3. The set S d is bounded i.e. ∀d > 0, ∀vw ∈ S d : Proof. For vw ∈ S d , by (3.6), Proof. By (3.6), The same can be shown for w by symmetry.
5 By Wedin's theorem (cf. matrix perturbation theory, see e.g. [28]), singular subspaces of M d associated with a positive singular value are continuously deformed with respect to d.

Multiplicative Updates for Nonnegative Factorization
In this section, the MU of Lee and Seung presented in Section 2.2 to find approximate solutions of (NMF) are generalized to (NF). Other than providing a way of computing approximate solutions of (NF), this result will also help us to understand why the updates of Lee and Seung are not very efficient in practice. The Karush-Kuhn-Tucker optimality conditions of the (NF) problem are the same as for (NMF) (see Section 2.2). Of course, any real matrix M can be written as the difference of two nonnegative matrices: M = P − N with P, N ≥ 0. This can be used to generalize the algorithm of Lee and Seung. In fact, (2.5) and (2.6) become and using the same idea as in Section 2.2, we get the following multiplicative update rules : Theorem 8. For V, W ≥ 0 and M = P − N with P, N ≥ 0, the cost function ||M − V W || F is nonincreasing under the following update rules: Proof. We only treat the proof for V since the problem is perfectly symmetric. The cost function can be split into m independent components related to each row of the error matrix, each depending on a specific row of P , N and V , which we call respectively p, n and v. Hence, we can treat each row of V separately, and we only have to show that the function is nonincreasing under the following update , ∀v 0 > 0. (3.16) F is a quadratic function so that . G has the following nice properties (see below): (1) G is an upper approximation of F i.e. G(v) ≥ F (v), ∀v; (2) The global minimum of G(v) is nonnegative and given by (3.16).
Therefore, the global minimum of G, given by (3.16), provides a new iterate which guarantee the monotonicity of F . In fact, It remains to show that (1) and (2) hold.
The global minimum of G is given by (3.16): As with standard multiplicative updates, convergence can be guaranteed with a simple modification: for any (V, W ) ≥ . Moreover, every limit point of this algorithm is a stationary point of the optimization problem (2.9).
Proof. We use exactly the same notation as in the proof of Theorem 3.15, so that is the sum of r independent quadratic terms, each depending on a single entry of v. Therefore, , and the monotonicity is proved. Let (V ,W ) be a limit point of a sequence {(V k , W k )} generated by (3.17). The monotonicity implies that {||M − V k W k || F } converges to ||M −VW || F since the cost function is bounded from below.
In order to implement the updates (3.15), one has to choose the matrices P and N . It is clear that ∀P, N ≥ 0 such that M = P − N , there exists a matrix C ≥ 0 such that the two components P and N can be written P = M + + C and N = M − + C. When C goes to infinity, the above updates do not change the matrices V and W , which seems to indicate that smaller values of C are preferable. Indeed, in the case r = 1, one can prove that C = 0 is an optimal choice: Proof. The second inequality of (3.19) is a consequence of Theorem 3.15. For the first one, we treat the inequality separately for each entry of v i.e. we prove that Let define v * i as the optimal solution of the unconstrained problem i.e.
Moreover, 1 ≤ a+d e+b+d since v 2i is a better solution than v i (Theorem 3.15). Finally, Unfortunately, this result does not hold for r > 1. This is even true for nonnegative matrices, i.e. one can improve the effect of a standard Lee and Seung multiplicative update by using a well-chosen matrix C. However, in practice, it seems that the choice of a proper matrix C is nontrivial and cannot accelerate significantly the speed of convergence.

How good are the Multiplicative Updates (MU) of Lee and Seung?
In this section, we use Theorem 8 to interpret the multiplicative rules for (NMF) and show why the HALS algorithm performs much better in practice.

An Improved Version of the MU
The aim of the MU is to improve a current solution (V, W ) ≥ 0 by optimizing alternatively V (W fixed), and vice-versa. In order to prove the monotonicity of the MU, ||M − V W || F was shown to be nonincreasing under an update of a single row of V (resp. column of W ) since the objective function can be split into m (resp. n) independent quadratic terms, each depending on the entries of a row of V (resp. column of W ); cf. proof of Theorem 8. However, there is no guarantee, a priori, that the algorithm is also nonincreasing with respect to an individual update of a column of V (resp. row of W ). In fact, each entry of a column of V (resp. row of W ) depends on the other entries of the same row (resp. column) in the cost function. The next theorem states that this property actually holds.
, ∀k, (4.1) i.e. under the update of any column of V or any row of W using the MU (2.7).
Proof. This is a consequence of Theorem 8 using P = M and N = i =k V :i W i: .
Corollary 2 sheds light on a very interesting fact: the multiplicative updates are also trying to optimize alternatively the columns of V and the rows of W using a specific cyclic order: first, the columns of V and then the rows of W . We can now point out two ways of improving the MU: 1. When updating a column of V (resp. a row of W ), the columns (resp. rows) already updated are not taken into account: the algorithm uses their old values; 2. The multiplicative updates are not optimal: P = M + and N = M − (cf. Theorem 10). Moreover, there is actually a closed-form solution for these subproblems (cf. HALS algorithm, Section 2.3).
Therefore, using Theorem 10, we have the following new improved updates

2)
with R k = M − i =k V :i W i: . Moreover, the updates (4.2) perform better than the updates (4.1), but worse than the updates (2.12-2.13) which are optimal.
Proof. This is a consequence of Theorem 8 and 10 using P = (R k ) + and N = (R k ) − . In fact, ||M − V W || F = ||R k − V :k W k: || F . Figure 3 shows an example of the behavior of the different algorithms: the original MU (Section 2.2), the improved version (Corollary 3) and the optimal HALS method (Section 2.3). The test was carried out on a commonly used data set for NMF: the cbcl face database 6 ; 2429 faces (columns) consisting each of 19 × 19 pixels (rows) for which we set r = 40 and we used the same scaled (see Remark 3 below) random initialization and the same cyclic order (same as the MU i.e. first the columns of V then the rows of W ) for the three algorithms. We observe that the MU converges significantly less Figure 3: Comparison of the MU of Lee and Seung (2.7), the optimal rank-one NF multiplicative updates (4.2) and the HALS (2.12-2.13) applied to the cbcl face database.

An example
rapidly than the two other algorithms. There do not seem to be good reasons to use either the MU or the method of Corollary 3 since there is a closed-form solution (2.12-2.13) for the corresponding subproblems. Finally, the HALS algorithm has the same computational complexity [16] and performs provably much better than the popular multiplicative updates of Lee and Seung. Of course, because of the NP-hardness of (NMF) and the existence of numerous locally optimal solutions, it is not possible to give a theoretical guarantee that HALS will converge to a better solution than the MU: although its iterations are locally more efficient, they could still end up at a worse local optimum.
Remark 2. For r = 1, one can check that the three algorithms above are equivalent for (NMF). Moreover, they correspond to the power method [13] which converges to the optimal rank-one solution, given that it is initialized with a vector which is not perpendicular to the singular vector corresponding to the maximum singular value.
Remark 3. We say that (V, W ) is scaled if the optimal solution to the problem is equal to 1. Obviously, any stationary point is scaled ; the next Theorem is an extension of a result of Ho et al. [18].
Theorem 11. The following statements are equivalent (1) (V, W ) is scaled; (2) V W is on the boundary of B M 2 , 1 2 ||M || F , the ball centered at M 2 of radius 1 2 ||M || F ; The solution of (4.3) can be written in the following closed form is the scalar product associated with the Frobenius norm.
Since α = 1, so that (1) and (2) are equivalent. For the equivalence of (1) and (3), we have Theorem 11 can be used as follows: when you compute the error of the current solution, you can scale it without further computational cost. In fact, Note that the third term of (4.5) can be computed in O(max(m, n)r 2 ) operations since where ||A|| 1 = ij |A ij |. This is especially interesting for sparse matrices since only a small number of the entries of V W (which could be dense) need to be computed to evaluate the second term of (4.5).

Biclique Finding Algorithm
In this section, an algorithm for the maximum edge biclique problem whose main iteration requires O(|E|) operations is presented. It is based on the multiplicative updates for Nonnegative Factorization and the strong relation between these two problems (Theorems 4, 6 and 7). We compare the results with other algorithms with iterates requiring O(|E|) operations using the DIMACS database and random graphs.

Description
For d sufficiently large, stationary points of (NF-1d) are close to bicliques of (MBP) (Theorem 7). Moreover, the two problems have the same cost function. One could then think of applying an algorithm that finds stationary points of (NF-1d) in order to localize a large biclique of the graph generated by M b . This is the idea of Algorithm 1 using the multiplicative updates (3.15) with A priori, it is not clear what value d should take. Following the spirit of homotopy methods, we chose to start the algorithm with a small value of d and then to increase it until the algorithm converges to a biclique of (MBP).
We observed that initial value of d should not be chosen too large: otherwise, the algorithm often converges to the trivial solution: the empty biclique. In fact, in that case, the denominators in (5.1) and (5.2) will be large, even during the initial steps of the algorithm, and the solution is then forced to converge to zero. Moreover, since the denominators in (5.1) and (5.2) depend on the graph density, the denser the graph is, the greater d 0 can be chosen and vice versa. On the other hand, since our algorithm is equivalent to the power method for d = 0 (cf. Remark 2), if d 0 is chosen too small, it will converge to the same solution: the one initialized with the best rank-one approximation of M b .
For the stopping criterion, one could, for example, wait until the rounding of vw coincides with a feasible solution of (MBP).

Other Algorithms in O(|E|) operations
We briefly present here two other algorithms to find maximal bicliques using O(|E|) operations per iteration.

Motzkin-Strauss Formalism
In [12], the generalized Motzkin-Strauss formalism for cliques is extended to bicliques by defining the optimization problem max x α i = 1}, F β y = {y ∈ R n + | n i=1 y β i = 1} and α, β 1. Nonincreasing multiplicative updates for this problem are then provided: This algorithm does not necessarily converge to a biclique: if α and β are not sufficiently small, it may converge to a dense bipartite subgraph (a bicluster). In fact, for α = β = 2, it converges to an optimal rank-one solution of the unconstrained problem as our algorithm does for d = 0. In [12], it is suggested to use α and β around 1.05. Finally, α = β will favor one side of the biclique. We will use α = β.

Greedy Heuristic
The simplest heuristic one can imagine is to add, at each step, a vertex which is connected to most vertices in the other side of the bipartite graph. Each time a vertex is selected, the next choices are restricted in order to get a biclique eventually: the vertices which are not connected to the one you have just chosen are deleted. The procedure is repeated on the remaining graph until you get a biclique. One can check that this produces a maximal biclique.

Results
We first present some results for graphs from the DIMACS graph dataset 7 . We extracted bicliques in those (not bipartite) graphs using the preceding algorithms. We performed 100 runs, 200 iterations each, for the two algorithms with the same initializations. We tried to choose appropriate parameters 8 for both algorithms. Table 1 displays the cardinality of the biclique extracted by the different algorithms. Table 2 shows the results for random graphs: we have generated randomly 100 graphs with 100 vertices for different densities (the probability of an edge to belong to the graph is equal to the density). The average numbers of edges in the solutions for the different algorithms are displayed for each density. We kept the same configuration as for the DIMACS graphs (same initializations, 100 runs for each graph, 200 iterations). It seems that the multiplicative updates generates, in general, better solutions, especially when dealing with dense graphs. The algorithm based on the Motzkin-Strauss formalism seems less efficient and is more sensitive to the choice of its parameters.   • It is applicable to non-binary matrices i.e. weighted graphs.
• It is possible to favor one side of the biclique. In fact, the multiplicative updates for NF can be adapted using the same developments as in Section 3.3 to cost functions with regularization terms, e.g. min v,w≥0 ||M − vw|| 2 F + α||v|| 2 2 + β||w|| 2 2 .
• If d is kept sufficiently small, for example replacing d = αd by d = min(αd, d m ) for some d m > 0, there is no guarantee that the algorithm will converge to a biclique. However, the negative entries in M d will enforce the corresponding entries of the solutions of (NF-1d) to be small (recall that Theorem 7 states that, for d sufficiently large, they will be equal to zero). Therefore, by rounding these solutions, instead of a biclique, one gets a dense submatrix of M b i.e. a bicluster. Algorithm 1 can then be used as a biclustering algorithm. The density of the corresponding submatrix will depend on the choice of d m . Table 3 gives an example of such behavior.

Conclusion
We have introduced Nonnegative Factorization (NF), a new variant of Nonnegative Matrix Factorization (NMF), and proved its NP-hardness for any fixed rank by reduction to the maximum edge biclique problem. The multiplicative updates for NMF can be generalized to NF and provide a new interpretation of the algorithm of Lee and Seung, which explains why it does not perform well in practice. We also developed an heuristic algorithm for the biclique problem whose iterations require O(|E|) operations, based on theoretical results about stationary points of a specific rank-one nonnegative factorization problem (NF-1d) and the use of multiplicative updates.
To conclude, we point out that none of the algorithms presented in this paper is guaranteed to converge to a globally optimal solution (and, to the best of our knowledge, such an algorithm has not been proposed yet) ; this is in all likelihood due to the NP-hardness of the NMF and NF problems. Indeed, only convergence to a stationary point has been proved for the algorithms of Sections 2 and 3, a property which, while desirable, provides no guarantee about the quality of the solution obtained (for example, nothing prevents these methods from converging to a stationary but rank-deficient solution, which in most cases could be further improved). Finally, no convergence proof for the biclique finding algorithm introduced in Section 5 is provided (convergence results from the preceding sections no longer hold because of the dynamic updates of parameter d) ; however, this heuristic seems to give very satisfactory results in practice.