Non-Negative Matrix Factorization Test Cases

I Introduction

^†^†footnotetext: This material is based in part upon work supported by the NSF under grant number DMS-1312831. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

What do document clustering, recommender systems, and audio signal processing have in common? All of them are problems that involve finding patterns buried in noisy data. As a result, these three problems are common applications of algorithms that solve non-negative matrix factorization, or NMF [2, 6, 7].

Non-negative matrix factorization involves factoring some matrix $A$ , usually large and sparse, into two factors $W$ and $H$ , usually of low rank

A = W H

(1)

Because all of the entries in $A$ , $W$ , and $H$ must be non-negative, and because of the imposition of low rank on $W$ and $H$ , an exact factorization rarely exists. Thus NMF algorithms often seek an approximate factorization, where $W H$ is close to $A$ . Despite the imprecision, however, the low rank of $W$ and $H$ forces the solution to describe $A$ using fewer parameters, which tends to find underlying patterns in $A$ . These underlying patterns are what make NMF of interest to a wide range of applications.

In the decades since NMF was introduced by Seung and Lee [5], a variety of algorithms have been published that compute NMF [1]. However, the non-deterministic nature of these NMF algorithms make them difficult to test. First, NMF asks for approximations rather than exact solutions, so whether or not an output is correct is somewhat subjective. Although cost functions can quantitatively indicate how close a given solution is to being optimal, most algorithms do not claim to find the globally optimal solution, so whether or not an algorithm gives useful solutions can be ambiguous. Secondly, all of the algorithms produced so far are stochastic algorithms, so running the algorithm on the same input multiple times can give different outputs if they use different random number sequences. Thirdly, the algorithms themselves, though often simple to implement, can have very complex behavior that is difficult to understand. As a result, it can be hard to determine whether a proposed algorithm really “solves” NMF.

This paper proposes some test cases that NMF algorithms should solve verifiably. The approach uses very simple input, such as matrices that have exact non-negative factorizations, that reduce the space of possible solutions and ensure that the algorithm finds correct patterns with little noise. In addition, small perturbations of these simple matrices are also used, to ensure that small variations in the matrix $A$ do not drastically change the generated solution.

Ii Perturbations of order $ϵ$

Suppose NMF is applied to a non-negative matrix $A$ to get non-negative matrices $W$ and $H$ such that $A \approx W H$ . If $A$ is chosen to have an exact non-negative factorization, then the optimal solution satisfies $A = W H$ . Furthermore, if $A$ is simple enough, most “good” NMF algorithms will find the exact solution.

For example, suppose $A_{0}$ is a non-negative square diagonal matrix, and the output $W_{0}$ and $H_{0}$ is also specified to be square. Let the diagonal $n \times n$ matrix $A_{0}$ be denoted $A_{0} = diag (a_{0})$ , where $a_{0}$ is an $n$

-dimensional vector, so that the diagonal entries

A_{0} (i, i)

are

a_{0} (i)

. It is easy to show that

W_{0}

and

H_{0}

must be monomial matrices (diagonal matrices under a permutation) [3]. Ignoring the permutation and similarly denoting

W_{0} = diag (w_{0})

and

H_{0} = diag (h_{0})

, then

a_{0} (i) = w_{0} (i) h_{0} (i)

for applicable

i

. Such diagonal matrices

A_{0}

were given as input to the known NMF algorithms described in the next section, and all of the algorithms successfully found exact solutions in the form of monomial matrices for

W_{0}

and

H_{0}

One way to analyze the properties of an algorithm is to perturb the input by a small amount $ϵ > 0$ and see how the output changes. Formally, if the input $A_{0}$ gives output $W_{0} H_{0}$ , then the output generated from $A_{0} + ϵ A_{1}$ can be approximated as $(W_{0} + ϵ W_{1}) (H_{0} + ϵ H_{1})$ . It is assumed that $ϵ$ is sufficiently small that $ϵ^{2}$ terms are negligible.

For the test case, the nonzero entries of $A_{1}$ were chosen to be the on the superdiagonal (the first diagonal directly above the main diagonal). This matrix is denoted as $A_{1} = diag (a_{1}, 1)$ , where $a_{1}$ is an $n - 1$ -dimensional vector such that $A_{1} (i, i + 1) = a_{1} (i)$ . The resulting matrix $A_{0} + ϵ A_{1}$ has $O (1)$ entries on its main diagonal, $O (ϵ)$ entries on the superdiagonal, and zeroes elsewhere. It is assumed that all the vector entries $a_{0} (i)$ and $a_{1} (i)$ are of comparable magnitude.

Iii Results from Various Algorithms

Three published NMF algorithms were implemented and run with input of the form $A = A_{0} + ϵ A_{1}$ as described above. Algorithm 1 was the multiplicative update algorithm described by Seung and Lee in their groundbreaking paper [5], which was run for $10^{6}$ iterations in each test. Algorithm 2 was the ALS algorithm described in [1], and which was run for $10^{6}$ iterations as well. Algorithm 3 was a gradient descent method as described by Guan and Tao [4], and was run for $10^{4}$ iterations. These three algorithms were chosen because they were representative and easy-to-implement algorithms of three distinct types. Many published NMF algorithms are variations of these three algorithms.

The experiments began with the simplest nontrivial case, in which $A$ is a $2 \times 2$ matrix with only three nonzero entries, with fixed $a_{0} = [11]$ and $a_{1} = [1]$ , while $ϵ$ was varied over several different values. Each of the algorithms used randomness in the form of initial seed values for $W$ and $H$ . The random seeds were held constant as $ϵ$ varied. As a result, the outputs from the algorithms with different values of $ϵ$ were comparable within each test case.

For the $2 \times 2$ case, it is possible to enumerate all of the non-negative exact factorizations of $A$ . Given that the factors $W$ and $H$ are also $2 \times 2$ matrices, they can be written as shown below.

[\begin{matrix} m & n p & q \end{matrix}] [\begin{matrix} r & s t & u \end{matrix}] = [\begin{matrix} 1 & ϵ 1 \end{matrix}]

(2)

Multiplying the matrices directly produces the the following four equations:

$m r + n t$	$= 1$	(3)
$m s + n u$	$= ϵ$	(4)
$p r + q t$	$= 0$	(5)
$p s + q u$	$= 1$	(6)

Recall that all entries must be non-negative, so from equation (5), either $p$ or $r$ must be 0, and either $q$ or $t$ must be 0. Furthermore, it cannot be that $p = q = 0$ because that would contradict equation (6), and it cannot be that $r = t = 0$ because that would contradict equation (3). Thus two cases remain: $p = t = 0$ and $q = r = 0$ .

Substituting $p = t = 0$ into equations (3), (4), and (6) and solving for $r$ , $s$ , and $u$ gives

r = \frac{1}{m}, s = \frac{1}{m} (ϵ - \frac{n}{q}), u = \frac{1}{q}

(7)

Likewise, substituting $q = r = 0$ into (3), (4), and (6) and solving for $s$ , $t$ , and $u$ to gives

(8)

Observe that these two solutions look similar. In fact, they differ merely by a permutation. In the first case, $W$ and $H$ have the same main diagonal and superdiagonal format as $A$ , and can be written in matrix notation as

W H = [\begin{matrix} w_{0} (1) & w_{1} (1) w_{0} (2) \end{matrix}] ⎡ ⎢ ⎣ \begin{matrix} \frac{1}{w_{0} (1)} & \frac{1}{w_{0} (1)} (ϵ - \frac{w_{1} (1)}{w_{0} (2)}) \frac{1}{w_{0} (2)} \end{matrix} ⎤ ⎥ ⎦

(9)

The second case can be written as $(W P) (P^{- 1} H)$ , where $P$ is the permutation matrix $[\begin{matrix} 11 \end{matrix}]$ .

All three of the algorithms tested gave solutions of this form 1000 times out of 1000, for each of several values of $ϵ$ . The consistency of the solutions enabled further analysis. The change in the solution can be measured by the change in the three parameters $w_{0} (1)$ , $w_{0} (2)$ , and $w_{1} (1)$ (ignoring the permutation if present). Figure 1 shows the change in each of the three parameters from the base case $A_{0}$ for several different values of $ϵ$ when input into Algorithm 1. Each of the values is the arithmetic mean of the corresponding values generated from 1000 different random seeds.

Fig. 1: The figure shows the slope associated with the change in each of the three parameters for each of several values of $ϵ$ . As $ϵ$ approaches zero on the right of the graph, the values of the slopes converge, showing that for sufficiently small $ϵ$ , each of the parameters is linear in $ϵ$ .

Of course, the precise values depend on the distribution of randomness used. But notice that as $ϵ$ approaches 0, the values of the three parameters become very nearly linear in $ϵ$ . The results for Algorithms 2 and 3 were very similar - they also showed linearity of the parameters in $ϵ$ , with comparable slopes.

However, $w_{1} (1)$ was not always linear in $ϵ$ , even for small $ϵ$ . In some cases, the difference approached 0 much more quickly. To see why this occurred, consider that the entries in $H$ could have been chosen to be the parameters rather than the entries in $W$ . Also, recall that in the base case $A_{0}$ , in which $ϵ = 0$ , $w_{1} (1) = h_{1} (1) = 0$ since both entries are off the diagonal. Thus, when either is linear in $ϵ$ , they are of the form $ϵ x$ for some slope $x$ . Since the solution is exact, it can be deduced that

w_{0} (1) h_{1} (1) + w_{1} (1) h_{0} (2) = ϵ

(10)

Therefore, in the cases that $w_{1} (1)$ approaches 0 very quickly, since $w_{0} (1)$ approaches a large, stable value as $ϵ$ approaches 0, $h_{1} (1)$ must be nearly linear in $ϵ$ . So in the cases that $w_{1} (1)$ is not linear in $ϵ$ , its symmetrical counterpart, $h_{1} (1)$ , is. To simplify this complication out of the data, the parameters in $W$ were chosen when $w_{1} (1)$ was closer to linearity in $ϵ$ , and the parameters in $H$ were chosen when $h_{1} (1)$ was closer to linearity in $ϵ$ .

Curiously, although it was possible for $w_{1} (1)$ and $h_{1} (1)$ to “split” the nonlinearity so that both were somewhat linear, this rarely occurred. All three algorithms preferred to make one of them very close to linear at the expense of the other. When $w_{1} (1)$ approached zero very rapidly, by equations (3) and (4), $h_{1} (1) = ϵ h_{0} (1)$ , and similarly, when $h_{1} (1)$ is negligible, $w_{1} (1) = ϵ h_{0} (2)$ .

Next, different values for the entries of $a_{0}$ and $a_{1}$ were tried, so they had a range of entries rather than all 1’s. The algorithms all behaved similarly; up to permutation, they satisfied the following formula

W H = [\begin{matrix} w_{0} (1) & w_{1} (1) w_{0} (2) \end{matrix}] ⎡ ⎢ ⎢ ⎣ \begin{matrix} \frac{a_{0} (1)}{w_{0} (1)} & \frac{a_{1} (1)}{w_{0} (1)} (ϵ - \frac{w_{1} (1) a_{0} (2)}{a_{1} (1) w_{0} (2)}) \frac{a_{0} (2)}{w_{0} (2)} \end{matrix} ⎤ ⎥ ⎥ ⎦

(11)

Note that equation (9) is just a special case of this equation in which $a_{0} (1) = a_{0} (2) = a_{1} (1) = 1$ . The same phenomena was also observed in which the algorithm usually made one of $w_{1} (1)$ and $h_{1} (1)$ be nearly linear in $ϵ$ and the other approach zero rapidly, rather than having both entries be non-negligible. As long as the entries of $a_{0}$ and $a_{1}$ are roughly on the order of 1, the algorithms operated similarly.

The next case examined set $A$ to be a $3 \times 3$ matrix. Using similar logic to the $2 \times 2$ case, it can be deduced that any exact factorization of $A$ is likely to be of the form

⎡ ⎢ ⎣ \begin{matrix} w_{0} (1) & w_{1} (1) w_{0} (2) & w_{1} (2) w_{0} (3) \end{matrix} ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ \begin{matrix} h_{0} (1) & h_{1} (1) h_{0} (2) & h_{1} (2) h_{0} (3) \end{matrix} ⎤ ⎥ ⎦

(12)

Indeed, all three algorithms always gave solutions of this form. In fact, most of the time there were two more zero entries than necessary - either $w_{1} (1)$ or $h_{1} (1)$ , and either $w_{1} (2)$ or $h_{1} (2)$ . This is similar to the way that $w_{1} (1)$ or $h_{1} (1)$ often approached 0 rapidly in the $2 \times 2$ case. To note another similarity to the $2 \times 2$ case, whenever $w_{1} (i)$ was significant and $h_{1} (i)$ was not, $w_{1} (i)$ was very close to $ϵ w_{0} (i + 1)$ - in similar situations $h_{1} (i)$ was approximately $ϵ h_{0} (i)$ .

As a result, there were 4 distinct configurations of the nonzero elements in the solutions, as given by Figure 2.

Type	equal to 0
Type I	$w_{1} (1), w_{1} (2)$
Type II	$h_{1} (1), h_{1} (2)$
Type III	$w_{1} (1), h_{1} (2)$
Type IV	$h_{1} (1), w_{1} (2)$

Fig. 2: We categorized the solutions when

A

was a

3 \times 3

matrix by where the non-negligible entries in the solution were. For each type, this table shows which entries that are usually positive are negligible.

Note that Type IV appears to be an inexact solution; since it has positive $w_{1} (1)$ and $h_{1} (2)$ , the entry at position $A (1, 3) = w_{1} (1) h_{1} (2)$ in the product $W H$ would have to be nonzero. However, both $w_{1} (1)$ and $h_{1} (2)$ , like all entries on the superdiagonal, are $O (ϵ)$ , so $w_{1} (1) h_{1} (2)$ is $O (ϵ^{2})$ , and is considered negligible. In fact, most of the solutions generated by the algorithms had nonzero values for entries that were supposed to be zero, but for this analysis anything below $O (ϵ^{2})$ was considered negligible.

Each algorithm was run 100 times on the $3 \times 3$ input with $w_{0} = [111]$ , $w_{1} = [11]$ , and $ϵ = 10^{- 3}$ . The solutions were categorized by the solution type in Figure 2. The distributions of the solutions by algorithm type are given in Figure 2. Note that some solutions did not have two negligible entries among $w_{1} (1)$ , $w_{1} (2)$ , $h_{1} (1)$ , and $h_{1} (2)$ , in which case the smaller entry was ignored for the sake of sorting - this accounted for about 20% of the three algorithms, the majority occurring in Algorithm 1.

Fig. 3: Categorized the solutions for $A$ being a $3 \times 3$ matrix by where the non-negligible entries in the solution were. This chart shows how often each algorithm generated a solution of each type out of 100 cases. Type II (in which $H$ is diagonal) was the most common among all the algorithms, but by differing amounts.

It is significant to note that even the solutions that didn’t fall cleanly into a “type” still satisfied the pattern given in (12). It seems that an NMF algorithm should satisfy this pattern, but little more is required.

Next, entries in $a_{0}$ and $a_{1}$ , were changed as in the $2 \times 2$ case. As long as the entries were $O (1)$ (as opposed to $O (ϵ)$ or $O (\frac{1}{ϵ})$ ), the behavior of the algorithms was similar.

Finally, $A$ larger than $3 \times 3$ were examined. Several different sizes of matrices were tested, ranging from $4 \times 4$ to $20 \times 20$ , always keeping $A$ , $W$ , and $H$ square, with positive entries only on the main diagonal and the superdiagonal. The experiments followed the same general pattern; nonzero entries in $W$ and $H$ appeared only on the main diagonal and superdiagonal. Using similar logic to the $2 \times 2$ and $3 \times 3$ cases, it can be shown that these are the only exact solutions. However, in practice, as the matrices get larger, exceptions to this pattern become more common, particularly in Algorithm 3. The general rule seems to mostly hold (over half the time) until $A$ becomes around $20 \times 20$ . Note, however, that because the run-time of the algorithms are cubic in the size of the matrix, at best, the sample size for large matrices is small.

Iv Proposed Tests for NMF Algorithms

Since all three algorithms, which cover a variety of approaches to NMF, had a lot in common in their solutions, it is propose that these inputs $A$ could be used as a test case of an NMF algorithm implementation. In this section, it is proposed how such test cases could be executed.

The test begins with input of the form

A = A_{0} + ϵ A_{1} = diag (a_{0}) + ϵ diag (a_{1}, 1)

(13)

$A$ is square, and preferably somewhere between $3 \times 3$ and $8 \times 8$ in size, although bigger inputs may be useful as well. The entries should vary between tests. Each test should start by using $ϵ = 0$ so that $A$ is diagonal. The results of this test should have $W$ and $H$ monomial - only one nonzero element in each row and column. Ignore entries that are below $O (10^{- 10})$ , for the entirety of testing, as any such entries are negligible.

If $W$ or $H$ is not monomial, or if the product $W H$ is not equal to $A$ to within a negligible margin of error, the algorithm fails the test. Otherwise, the generated solution can be used to find the permutation matrix $P$ that makes $W P$ and $P^{- 1} H$ diagonal by replacing the nonzero entries of $H$ with 1’s. Since $A = W H$ is diagonal, $W P$ is also diagonal, and since $I = P^{- 1} P$ is diagonal, so is $P^{- 1} H$ . Knowing $P$ will make the rest of the testing much simpler since it is easier to identify whether a solution is of the form given above when it is not permuted.

Next, run the test again using a positive value for $ϵ$ ; $ϵ = 10^{- 3}$ seems to work well, although using a variety of $ϵ$ is also recommended. Make sure to use the same random seeds that were used in the $ϵ = 0$ test to produce corresponding output. Then check that the $W$ and $H$ given by the algorithm are such that $W P$ and $P^{- 1} H$ have nonzero entries only on the two diagonals that they are supposed to. If this doesn’t hold, changing $ϵ$ might have changed which permutation returns $W$ and $H$ to the proper form, so check again; this happens more commonly among larger matrices than smaller ones. However, if $W$ and $H$ really do break the form, or $A \neq W H$ , the algorithm fails the test on this input. Otherwise, it passes.

Note that even widely accepted algorithms do fail these tests occasionally, especially with matrices larger than $8 \times 8$ , so it’s advisable to perform the test many times to get a more accurate idea of an algorithm’s performance.

V Conclusion

This paper proposes an approach to the problem of testing NMF algorithms by running the algorithms on simple input that can produce an exact non-negative factorization, and perturbations of such input. In particular, square matrices with $O (1)$ entries on the main diagonal and $O (ϵ)$ entries on the superdiagonal are proposed, because they have exact solutions that can enumerated mathematically, or because they are perturbations of matrices with exact solutions.

The test cases have been used as input on three known NMF algorithms that represent a variety of algorithms, and all of them behaved similarly, which suggests testable, quantifiable behaviors that many NMF algorithms share. These test cases offer one approach for testing candidate NMF implentations to help determine whether it behaves as it should.

Acknowledgment

The authors would like to thank Dr. Alan Edelman for providing and overseeing this research opportunity, and Dr. Vijay Gadepally for his advice and expertise.

References

[1] Berry Browne Langville Pauca and Plemmons, Algorithms and applications for approximate nonnegative matrix factorization, Computational Statistics and Data Analysis 52 (2007), 155–173.
[2] Rainer Gemulla, Erik Nijkamp, Peter J. Haas, and Yannis Sismanis, Large-scale matrix factorization with distributed stochastic gradient descent
, Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (New York, NY, USA), KDD ’11, ACM, 2011, pp. 69–77.
[3] John Gilbert, personal communication, Sep 2015.
[4] N. Guan, D. Tao, Z. Luo, and B. Yuan, Nenmf: An optimal gradient method for nonnegative matrix factorization, IEEE Transactions on Signal Processing 60 (2012), no. 6, 2882–2898.
[5] Daniel D. Lee and H. Sebastian Seung, Algorithms for non-negative matrix factorization, Advances in Neural Information Processing Systems 13 (T. K. Leen, T. G. Dietterich, and V. Tresp, eds.), MIT Press, 2001, pp. 556–562.
[6] Suvrit Sra and Inderjit S. Dhillon, Generalized nonnegative matrix approximations with bregman divergences, Advances in Neural Information Processing Systems 18 (Y. Weiss, B. Schölkopf, and J. C. Platt, eds.), MIT Press, 2006, pp. 283–290.
[7] Wenwu Wang, Instantaneous vs. convolutive non-negative matrix factorization: Models, algorithms and applications, Machine Audition: Principles, Algorithms and Systems: Principles, Algorithms and Systems (2010), 353.

Non-Negative Matrix Factorization Test Cases

Authors

Non-negative matrix factorization with sparseness constraints

Fast Conical Hull Algorithms for Near-separable Non-negative Matrix Factorization

PL-NMF: Parallel Locality-Optimized Non-negative Matrix Factorization

Integrating Hypertension Phenotype and Genotype with Hybrid Non-negative Matrix Factorization

Semi-Orthogonal Non-Negative Matrix Factorization

Relative Pairwise Relationship Constrained Non-negative Matrix Factorisation

A sampling algorithm to compute the set of feasible solutions for non-negative matrix factorization with an arbitrary rank

I Introduction

Ii Perturbations of order $ϵ$

Iii Results from Various Algorithms

Iv Proposed Tests for NMF Algorithms

V Conclusion

Acknowledgment

References

Non-Negative Matrix Factorization Test Cases

Authors

Related Research

Non-negative matrix factorization with sparseness constraints

Fast Conical Hull Algorithms for Near-separable Non-negative Matrix Factorization

PL-NMF: Parallel Locality-Optimized Non-negative Matrix Factorization

Integrating Hypertension Phenotype and Genotype with Hybrid Non-negative Matrix Factorization

Semi-Orthogonal Non-Negative Matrix Factorization

Relative Pairwise Relationship Constrained Non-negative Matrix Factorisation

A sampling algorithm to compute the set of feasible solutions for non-negative matrix factorization with an arbitrary rank

This week in AI

I Introduction

Ii Perturbations of order ϵ

Iii Results from Various Algorithms

Iv Proposed Tests for NMF Algorithms

V Conclusion

Acknowledgment

References

Ii Perturbations of order $ϵ$