Non-Negative Matrix Factorization with Sinkhorn Distance

1 Introduction

Data representation plays a fundamental role in various problems in the fields of pattern recognition, information retrieval and computer vision [Bengio et al., 2013]. A good representation can significantly improve the performance of algorithms. Nowadays, it is easy to collect massive data with lots of features for the problems, which leads to large-scale high-dimensional data set. However, the intrinsic degrees of freedom could be far less. Naturally, one might hope to find the space where the data can be represented by a small number of semantic concepts. To achieve this goal, a large number of methods were proposed in the last decades, among which matrix factorization based methods [Wall et al., 2003; Srebro et al., 2004] have received considerable attentions.

Since its psychological and physiological interpretation of naturally occurring data, Non-negative Matrix Factorization (NMF) [Lee and Seung, 1999] is especially striking among these methods. Different from other methods, NMF uses nonnegativity constraints and these constraints lead to a partsbased representation, which may be as same as the form that data is represented in human brains[Lee and Seung, 1999]. It has been shown to be superior to other matrix factorization based methods in various applications, such as face recognition [Guan et al., 2012] and document clustering [Shahnaz et al., 2006; Xu et al., 2003]. One shortcoming of the original NMF is that it ignores the intrinsic structure of the data set. The problem is twofold. First, samples might be on a manifold and thus one may hope this information can be exploited to improve NMF's performance[Cai et al., 2011]. Second, features might correlate with each other, thus conventional L2 distance can not well measure the distance between samples [Sandler and Lindenbaum, 2011].

To address this limitation, researchers propose to improve NMF in different ways. [Cai et al., 2011] uses a nearest neighbor graph to model the local geometry structure and aims at finding a matrix factorization that preserves the graph structure. Therefore, the proposed GNMF keeps the local structure on the data manifold. Instead of constructing a local connection graph, EMD NMF [Sandler and Lindenbaum, 2011] uses Earth Mover's Distance (EMD) [Rubner et al., 2000] (also known as Wasserstein distance) to measure the difference between original data matrix and the product matrix. EMD is defined as the cost of the optimal transport plan for moving the mass between two histograms. It has been applied to a wide range of problems, including label propagation [Solomon et al., 2014], and supervised learning [Frogner et al., 2015]. EMD is sensitive to relationships between the different dimensions and agrees with perceptual dissimilarity better than other measures[Rubner et al., 2000]. Thus EMD NMF achieves more robust results over traditional L2 NMF for problems where the error mechanism follows complex local deformation[Sandler and Lindenbaum, 2011].

However, the computation of EMD NMF is time consuming, thus limits its scalability. Even though the wavelet EMD (WEMD)[Shirdhonkar and Jacobs, 2008], which is an accelerated approximation of EMD, is employed to accelerate optimization, EMD NMF is still slow. For a data matrix with 200 samples and 1024 features, one full iteration of WEMD based EMD NMF takes around 23 minutes and the algorithm converges after 4 hours by using Matlab on an Intel Core 2 Quad 2.5GHz processor[Shirdhonkar and Jacobs, 2008]. Moreover, for the WEMD-based EMD NMF algorithm, the factorization results may not be non-negative since only soft penalty is used to encourage non-negativity.

In this paper, we propose a novel NMF algorithm which exploits information in both data manifold and features correlation. Inspired by the recent progress on EMD [Cuturi, 2013; Frogner et al., 2015], we propose to model the relationship between different feature dimensions by unnormalized Wasserstein distance with entropic regularization (also named as Sinkhorn distance in short), whose gradient can be solved efficiently. Meanwhile, graph regularization based on this distance is incorporated to preserve local geometry structure. Our model can capture both manifold structure and features correlation and is named as Non-negative Matrix Factorization with Sinkhorn distance (SDNMF). We further propose an efficient multiplicative update algorithm for the proposed NMF model. The non-negative factorization can then be efficiently computed. Experimental results on real-world data sets demonstrate the effectiveness of our model and the efficiency of our proposed multiplicative update algorithm.

The rest of the paper is organized as follows: section 2 provides a brief review of EMD and its entropic regularization. Our SDNMF model and multiplicative update rules are introduced in section 3. Experiments on two real-world data sets are presented in section 4. Related work is presented in section 5. Finally, we provide some concluding remarks in section 6.

2 Preliminary

Notation. Here we briefly introduce the notations used in this paper. We use italic uppercase letters to denote matrices, bold lowercase letters to denote vectors. Given an arbitrary m⇥n matrix A,we define aj and aij as the j-th column vector and the (i↪ j)-th entry of matrix A respectively. ⊙ represents element-wise multiplication and ↵ represents element-wise division.

2.1 EarthMover's Distance

Given two normalized histograms x↪ y 2 Rm(Pt xt = Pt yt = 1), and a distance metric matrix M, the Earth Mover's Distance dM(x↪ y) [Rubner et al., 2000] is defined as:

Eq. (1) is the well-known transportation problem and the flow variable Tpq denotes the quantity transported from the p-th supply to the q-th demand. The parameter Mpq represents the ground distance between bins p and q. Usually Mpq is defined by L1 or L2 distance or is determined based on the priori knowledge of the features in the considered problem.

However, the computation of EMD is time consuming. [Cuturi, 2013] proposes a method to accelerate its computation. They smooth the classical optimal transportation problem with an entropic regularization term and show that the resulting optimum is also a distance which can be computed through Sinkhorn-Knopps matrix scaling algorithm at a speed that is several orders of magnitude faster than that of transportation solvers. They call the new distance as Sinkhorn distance and the the specific form is:

Moreover, [Frogner et al., 2015] proposes a relaxation that extends the smoothed transport to unnormalized measures. They replace the equality constraints on the transport marginals in Eq. (2) with soft penalties with respect to KL divergence and obtain an unconstrained approximate transport problem. The distance becomes:

where H(T) = - Pp↪q Tpq logTpq is the entropy of T. It can approximate the exact EMD closely with λ large enough.

Eq. (2), Eq. (3) and their gradients can be efficiently solved by Sinkhorn-like iterations [Frogner et al., 2015].

+ γ(g KL(T 1♣♣x) + g where g KL(w♣♣z) = wT log(w↵z)-1T w+1T z is the generalized KL divergence between w and z. It degenerates to Eq. (2) with γ large enough when Pt xt = Pt yt = 1.

3 SDNMF

Given a data matrix X = [xij ] = [x1↪ ▷▷▷↪ xn] 2 Rm⇥n , where each column is a sample vector, NMF tries to find two nonnegative matrices U = [uik] = [u1↪ ▷▷▷↪ ut] 2 Rm⇥t and V = [vjk] = [v1↪ ▷▷▷↪ vt] 2 Rn⇥t such that Y = [yij] = [y1 ↪ ▷▷▷↪ yn ] is an approximation of origin data matrix X under the new data representation yj = Ptk=1 ukvjk. Previous works usually assume the error between X and Y is due to Gaussian error. However, this assumption neglects the relationship between different dimensions of features. Considering the case where we have two documents represented by Bag-of-Words which have the same topic "soccer", one is represented as ¶"soccer" : 100↪ "football" : 0↪ ▷▷▷♦ and the other is represented as ¶"soccer" : 0↪ "football" : 100↪ ▷▷▷♦. It does not affect their topics if we change the count of "soccer" and the count of "football" as long as we keep the count of "soccer" and "football". Noted that the topic vector "soccer" produced by NMF will probably have nonzero count on both of them. It causes that the new representation produced by NMF will have nonzero count on both of them and the error of NMF increases. Obviously, this kind of error can not be modeled by Gaussian error. Therefore, we adopt EMD as the metric. Optimizing with respect to the exact EMD is costly and it's difficult to ensure Pi yij = Pi xij in the whole iterative solving process. Eq. (3) describes a regularized approximation which can be efficiently computed even with unnormalized data. So we change the objective function as:

In addition to the correlation of features, there exists a geometry structure in the sample space. We hope that the information of the geometric structure of samples can be exploited for better discovery of the basis [u1↪ ▷▷▷↪ ut]. Inspired by [Cai et al., 2011], we also use the manifold assumption on the data. This assumption can be explained as that if two data points xj , xs are close in the intrinsic geometry of the data manifold, then vj and vs, the representations of these two points in the new basis, are also close to each other. Formally, it can be written as:

where W is a nearest neighbor graph on a scatter of data points. Different from [Cai et al., 2011], our nearest neighbor graph is constructed based on Sinkhorn distance.

Combining both of them, we get our final objective function:

where ⇠ is a regularization parameter that controls the tradeoff between the features correlation and the data manifold.

3.1 Multiplicative Update Rules

The objective function O of SDNMF in Eq. (6) is not jointly convex in U and V . Therefore it is impractical to expect an algorithm to seek the global minimum of O. Fortunately, the objective function O is convex in U and V separately. Similar to [Lee and Seung, 2001], we also adopt a two-stage multiplicative update rules which can keep non-negativity and find a local minimum:

where T ⇤it s is the (i↪ t)-entry of the optimal transportationmatrix between xs and ys and can be solved by a slight variation of algorithm 1 in [Frogner et al., 2015] efficiently. What's more, the following theorem guarantees that the update rules of U and V in Eq. (7) and (8) converge and the final solution will be a local optimum. A detailed proof can be found in the appendix.

Theorem 1. The objective function O in Eq. (6) is nonincreasing under the update rules in Eq. (7) and (8). The objective function is permanent under these updates if and only if U and V are at a stationary point.

4 Experiments

In this section, we test the performance of the proposed algorithm in the context of two challenging tasks. Previous studies show that NMF is very powerful on clustering, especially in the document clustering [Shahnaz et al., 2006; Xu et al., 2003] and image clustering tasks [Guan et al., 2012]. For Image clustering, previous works usually align images first. But in the real world, images might be too complex to be aligned perfectly. There may still exist local deformation after alignment. In the first experiment, we use small random translation to simulate local deformation and show that our method is more robust than previous methods. In the second experiment, we use NMF to estimate texture descriptors for texture mosaic images and show that our method has a strong ability of description for images with high local variability.

4.1 Image Clustering with Translation Noise

Data Set. The data set used in this section is COIL20 image library [Nene et al., 1996], which contains well aligned 32⇥32 gray scale images of 20 objects viewed from 72 varying angles. However images might be too complex to be aligned perfectly in the real world. There will still exist local deformation after alignment. We use small random translation to simulate local deformation and test our method on the new data set to show its robustness. Specifically, our data sets are generated by:

1. Resizing each image to (32-δ)⇥(32-δ) and placing it at the center of a 32 ⇥ 32 blank image, where δ is a small integer that controls the degree of random translation noise.
2. Generating a random integer vector (i↪ j), where i and j are sampled from a discrete uniform distribution on ¶-δ↪ ▷▷▷↪ δ♦ independently.
3. Translating the resized image with a vector (i↪ j).

Some examples are shown in Figure 1.

Compared Algorithms. We compare our method with the following four popular clustering algorithms on data sets generated with δ = 0↪ 1↪ 2↪ 3↪ 4 respectively. For each compared

19.55◦0.80 19.52◦1.06 21.31◦0.77 18.70◦0.88 23.15◦1.56 27.88◦1.24 26.48◦1.08 29.79◦1.12 24.75◦0.51 34.28◦1.61

method, several parameter configurations are tested and the best performance is reported.

• ÷ Canonical K-means clustering method (Kmeans in short).
• ÷ Non-negative Matrix Factorization based clustering (NMF in short) [Lee and Seung, 2001].
• ÷ Graph Regularized Non-negative Matrix Factorization (GNMF in short) with Frobenius norm formulation. Following [Cai et al., 2011], we use binary weighting scheme for constructing the 5-nearest neighbor graph and set the regularization parameter λ to 100.
• ÷ Non-negative Matrix Factorization with Earth Mover's Distance (EMD NMF in short). Since the WEMD algorithm in [Sandler and Lindenbaum, 2011] is too slow (the detailed performance will be discussed later) and can not ensure its result is non-negative, we use our efficient algorithm to solve it 1 . We set the λ and γ to 100 and 10 respectively and use the 2D distance of pixels' location of the image as the ground metric.
• ÷ Non-negative Matrix Factorization with Sinkhorn Distance (SDNMF in short). We use binary weighting scheme for constructing the 5-nearest neighbor graph for its simplicity. We set the λ , γ and ⇠ to 100,1 and 10 respectively and use the 2D distance of pixels' location of the image as the ground metric.

We evaluate the clustering performance by comparing the obtained label of each sample with the label provided by the data set. Two metrics, the accuracy (AC) and the normalized mutual information (NMI) are used to measure the clustering result. The detailed definitions of these two metrics can be found in [Cai et al., 2005].

Clustering Results. Table 1 shows the clustering results on the COIL2 data sets with different levels of random translation noise. For each given random size δ, 10 test runs were conducted on different randomly chosen clusters. The table reports the mean and the standard error of the performance.

It can be seen from Table 1 that our method has the best performance in most cases. When the size of random noise is increasing, our method has the most robust performance.

Recall that for a data matrix with 200 samples and 1024 features, one full iteration of the algorithm in [Sandler and Lindenbaum, 2011] takes around 23 minutes and the algorithm converges after 4 hours by using Matlab on an Intel Core 2 Quad 2.5 GHz processor. Our algorithm needs only

1 Setting ⇠ to zero, our objective function Eq. (6) is another alternative approximation of EMD NMF

2.9 seconds per full iteration and converges after 290.9 seconds by using Matlab on an Intel i7-4790K 4.0GHz processor on a data set with 1440 samples and 1024 features. It demonstrates the efficiency of our algorithm.

Parameters Selection. Our SDNMF model has three essential parameters: the regularization parameter of Sinkhron distance λ, the relaxation parameter of unnormalized Sinkhorn distance γ and the trade-off parameter ⇠. Figure 2 shows how the average performance of SDNMF varies with these parameters. As we can see, the performance of SDNMF is very stable with respect to these parameters. SDNMF achieves consistently good performance when these parameters varies in a large range.

4.2 Texture descriptor estimation

In this section we do the same experiment in section 6.2 of [Sandler and Lindenbaum, 2011]. Our task is estimating the texture descriptors associated with each texture class of the mosaic from [Haindl and Mikeˇs, 2008] (some examples are shown in Figure 3). Moreover, roughly classification of the textures in each mosaic location (e.g., for consecutive segmentation) is also expected. [Haindl and Mikeˇs, 2008] contains online generated 512⇥512 mosaics with different numbers of textures and the number can be chosen from 3 to 12. For each number, we choose 9 samples to make up our data set. For our task, we consider the texture in non-overlapping square image patches and assume the texture in each block is a positive mixture of the basic textures. Then NMF is used to analyse it.

For each texture, we assume there exists a vector descriptor utrue i connected with it. Therefore, a mosaic with C textures is connected with a texture descriptor matrix U true = [utrue 1 ↪ ▷▷▷↪ utrue C ]. We expect that the texture descriptor in the j-th image patch should be xj = U truevtrue j , where vtrue j is the vector of true fractions of the j-th block area associated with each texture class. Then we can use NMF to roughly estimate texture descriptors and classify the textures in each mosaic location.

The specific process is:

1. Converting the image to a new representation which each location is represented by a vector of Gabor responses, since Gabor filters have been widely used in texture analysis [Ramakrishnan et al., 2002].
2. Dividing the image into N non-overlapping patches and computing the mean feature vector xj for each patches. Now, X = [x1↪ ▷▷▷↪ xN] is a new representation of the image.
3. Using NMF to find the factorization X ⇡ UV T .

The factorization results U = [u1↪ ▷▷▷↪ uC] and V = [v1↪ ▷▷▷vC] are the approximated representative texture descriptors and the approximated fraction of each texture in each patch respectively.

We conduct the tests with different number of patches N = 16, 64, 256, 1024 and these patches tessellate the image. Thus, N determines the patch size 128⇥128↪ 64⇥64↪ 32⇥32, and 16 ⇥ 16 pixels respectively. The U true consists of the mean descriptors on the single texture segments of the mosaic and the V true consists of the fractions of each class in the patch. We use average correlation between the true representation and the approximated one to evaluate the performance. This metric can be written as:

Figure 4: The performance of texture descriptor estimation

Figure 4 shows the performance of texture descriptor estimation. We still use NMF, GNMF and EMDNMF as compared algorithms. It can be seen from Figure 4 that our method has the best performance in most cases. Since textures exhibit lots of spatial variation, EMD might be preferable over L2. Therefore, EMD based method will be more robust when the variation increasing (patch size decreasing). What's more, L2 based graph regularizer will damage the performance.

5 RelatedWork

In the past decade, based on the standard NMF, many variants have been proposed to find effective data representation for various problems. The family of NMF algorithms have been applied in many areas, such as face recognition[Guan et al., 2012], clustering[Xu et al., 2003], hyperspectral unmixing[Jia and Qian, 2009], etc., and achieve remarkable performance.

The works that are closely related to our method could be roughly categorized into two lines. First, various difference metrics are adopted to measure the dissimilarity between the original data matrix and the product. [Kompass, 2007] proposes a generalized divergence measure that interpolates between square loss and Kullback-Leibler divergence for non-negative matrix factorization. [F´evotte et al., 2009] uses Itakura-Saito divergence to incorporate Bayesian priors. However, none of these works model the features correlation explicitly and they may fail to capture the relationships of different dimensions. [Sandler and Lindenbaum, 2011] first proposes to use Earth Mover's Distance as difference metric, and employs wavelet EMD approximation to derive practical algorithm. Since EMD is sensitive to the feature correlation and agrees with the perceptual dissimilarity better than other measures[Rubner et al., 2000], EMD NMF achieves more robust results for problems where the error mechanism follows complex local deformation. However, the WEMD based EMD NMF is still computationally costly. In contrast, both the Sinkhorn Distance that is considered in our method and its gradient can be computed efficiently by Sinkhorn-like matrix scaling algorithm.

Second, considering the local geometry structure of the data manifold, [Cai et al., 2011] constructs a nearest neighbor graph and finds a matrix factorization that preserves this graph structure. The proposed GNMF can keep the local structure on the data manifold, thus achieves superior results on clustering. GNMF is further extended to problems like coclustering[Shang et al., 2012], by adding more constraints. A semi-supervised NMF method is also proposed in [Liu et al., 2012] to incorporate available labeled samples, which increases the discriminative ability of the obtained representation. However, the local graphs in these works are based on conventional L2 distance or Kullback-Leibler divergence, which can not reflect the feature correlation. Meanwhile, currently we only consider the case without supervision.

6 Conclusions and FutureWork

In this paper we propose a new NMF method, SDNMF, which exploits knowledge in both data manifold and features correlation. We use EMD to utilize information of feature correlation and use a graph regularizer to keep the local geometric structure of the data manifold. Although optimizing with respect to the exact EMD objective function is computationally costly, our approximate Sinkhorn distance objective function is efficiently computed. By this approximation, we propose a fast implementation of the proposed algorithm. SDNMF outperforms previous NMF based algorithms in the context of two challenging computer vision tasks.

Although our multiplicative update algorithm is efficient, SDNMF has to pay high computational cost in graph construction and subsequent matrix calculation, which limits its applicability to large-scale problems. Inspired by [Liu et al., 2010] an interesting direction for future work may be to use landmark points to speed up our algorithm.

Acknowledgments

This work was supported by the National Basic Research Program of China (973 Program) under Grant 2012CB316400, National Natural Science Foundation of China (Grant nos. 61233011, 61125203), and National Youth Top-notch Talent Support Program.

A (Proofs of Theorem 1)

To prove Theorem 1, we will imitate the proof process from [Lee and Seung, 2001] and use its Definition 1 and Lemma 1. Please see [Lee and Seung, 2001] for details.

Definition 1. G(h↪ h0 ) is an auxiliary function for F(h) if the conditions

are satisfied.

Lemma 1. If G is an auxiliary function of F, then F is nonincreasing under the update

For simplicity, we only prove that O is non-increasing under the update step in Eq. (8) (the other half can be proved in similar way). Fixed U, we rewrite the objective function O as follows

According to Definition 1 and Lemma 1, we need to find an auxiliary function G(V↪ V 0 ) for it.

Lemma 2. Setting ↵iki = uikv(q) jk Pk uikv(q) jk and ↵ = [↵1k1 ↪ ▷▷▷↪ ↵mkm ]T . Then function

is an auxiliary function for F(V ).

Proof. It is straightforward to verify that G(V↪ V ) = F(V ). Since dλ↪γ M is convex [Frogner et al., 2015], we use the convexity for i = 1↪ ▷▷▷↪ m one by one and obtain

Thus we find a proper function G(V↪ V (q) ).

Then we can prove Theorem 1:

Proof of Theorem 1. Note that the gradient of dλ↪γ M (y↪ x) with respect to y is: [Frogner et al., 2015]

By setting the gradient of G(V↪ V (q) ) with respect to V to zero, we can find V (q+1) .

According to Definition 1 and Lemma 1, Theorem 1 is proved