Fast Dictionary Learning with a Smoothed Wasserstein Loss

Antoine Rolet
antoine.rolet@iip.ist.i.kyoto-u.ac.jp
Graduate School of Informatics,
Kyoto University

Marco Cuturi
Gabriel Peyrémcuturi@i.kyoto-u.ac.jp
gabriel.peyre@ceremade.dauphine.fr
Graduate School of Informatics,
CNRS, CEREMADE,
Kyoto University
Université Paris Dauphine

Abstract
We consider in this paper the dictionary
learning problem when the observations arenormalized histograms of features. This
problem can be tackled using non-negative
matrix factorization approaches, using typically Euclidean or Kullback-Leibler fitting errors. Because these fitting errors are separable and treat each feature on equal footing,they are blind to any similarity the features
may share. We assume in this work that wehave prior knowledge on these features. To
leverage this side-information, we propose to
use the Wasserstein (a.k.a. earth mover’s oroptimal transport) distance as the fitting error between each original point and its reconstruction, and we propose scalable algorithms to to so. Our methods build upon
Fenchel duality and entropic regularizationof Wasserstein distances, which improves not
only speed but also computational stability.
We apply these techniques on face images andtext documents. We show in particular that
we can learn dictionaries (topics) for bagof-word representations of texts using words
that may not have appeared in the original
texts, or even words that come from a different language than that used in the texts.

1

Introduction

Consider a collection X = (x1 , . . . , xm ) of m vectorsof dimension n. Learning a dictionary for X can bestated informally as the goal of finding k dictionary
elements D = (d1 , . . . , dk ) of the same dimension n
To appear in Proceedings of the 19th International Conference on Artificial Intelligence and Statistics (AISTATS)
2016, Cadiz, Spain. JMLR: W&CP volume 41. Copyright2016 by the authors.

such that each Xi can be reconstructed using such a
dictionary, namely such that there exists a matrix of
mixture weights Λ = (λ1 , . . . , λm ) such that X ' DΛ.
When all elements of X are non-negative, and if it isdesirable that all elements of D and Λ are nonnegative
too, this problem becomes that of Non-negative Matrix
Factorization (NMF) [Paatero and Tapper, 1994]. Lee
and Seung [2001] proposed two algorithms for NMF,
with the aim of solving problems of the form:
min

m
X

D∈Rn×k,Λ∈Rk×m++i=1`(xi , Dλi ) + R(D, Λ),

where ` is either the Kullback-Leibler divergence or thesquared Euclidean distance and R a regularizer. Dictionary learning and NMF have been used for various
machine learning and signal processing tasks, including (but not limited to) semantic analysis [Hofmann,
1999, Lee and Seung, 1999], matrix completion [Zhang
et al., 2006] and sound denoising [Schmidt et al., 2007].
Our goal in this paper is to generalize these approaches
using a regularized Wasserstein (a.k.a optimal transport [Villani, 2009] or earth mover’s [Rubner et al.,1998]) distance as the data fitting term `. Such distances can leverage additional knowledge on the spaceof features using a metric between features called theground metric. Since the seminal work of Rubner
et al. [1998], several hundred papers have successfully used EMD in applications. Some recent works
have for instance illustrated its relevance for textclassification [Kusner et al., 2015], image segmentation [Rabin and Papadakis, 2015] and shape interpolation [Solomon et al., 2015].
We motivate the idea of using a Wasserstein fitting error with a toy example illustrated in Figure 1. In thisexample we try to learn dictionaries for histogram representations of i.i.d. samples from mixtures of Gaussians. We consider n = 100 distributions ρ1 , . . . , ρn ,each of which is a mixture of three univariate Gaussians of unit variance, with centers picked independently using N (−6, 2), N (0, 2) and N (6, 2) respectively. The relative weights of these Gaussians areFast Dictionary Learning with a Smoothed Wasserstein Losspicked uniformly on [0, 1] and subsequently normalized
to sum to 1 for each distribution. We consider then asample of m observations for each distribution ρi , andrepresent each sample as a histogram xi of n = 100
bins regularly spaced on the segment [−12, 12]. Herethe features are points on the quantization grid, andthe ground metric is simply the Euclidean distance between these points. Wasserstein NMF recovers components which are centered around −6, 0 and 6 andresemble Gaussian pdfs. Because it is blind to themetric structure of R, KL NMF fail to recover such
intuitive components.

Data samples
0.1
x1
x2
x3

0.08
0.06
0.04
0.02
0
-6

0

6

Wasserstein NMF
0.1

KL NMF

0.1
d1
d2
d3

0.08
0.06

0.06

0.04

0.04

0.02

0.02

0

d1
d2
d3

0.08

0
-6

0

6

-6

0

6

Figure 1: Dictionaries learned on mixtures of threeGaussians with a random shift. Separable distances ordivergences do not quantify this noise well because it is
not additive in the space of histograms. Top: examplesof data histograms. Bottom: dictionary learned with
Wasserstein (left) and Kullback-Leibler (right) NMF.

Related Work Sandler and Lindenbaum [2009]were the first to consider NMF problems using a
Wasserstein loss.
They noticed that minimizing
Wasserstein fitting errors requires solving an extremely
costly linear program at each iteration of their blockcoordinate iteration. Because of this, they settle instead for an approximation of the Wasserstein distanceproposed by Shirdhonkar and Jacobs [2008]. However
this approximation can only be used when the features
are in Rd , and its complexity is exponential in d, making it practical only when d ≤ 3. Zen et al. [2014]
also proposed a semi-supervised method to learn D, Λand a ground metric parameter. Their approach is toalternatively learn the ground metric as proposed pre-viously in [Cuturi and Avis, 2014] and perform NMF
by solving two very high dimensional linear programs.
They apply their algorithm to histograms of small dimension (n ≤ 16).
Our Contribution The algorithms we propose tosolve dictionary learning and NMF problems with a
Wasserstein loss scale to problems with far more observations and dimensions than previously considered in
the literature [Sandler and Lindenbaum, 2009, Shirdhonkar and Jacobs, 2008]. This is enabled by an
entropic regularization of optimal transport [Cuturi,
2013] which results in faster and more stable computations. In contrast to previously considered approaches,
our approach can be applied with any ground metric.
We introduce this regularization in Section 2, and follow in Section 3 with a detailed presentation of our algorithms for Wasserstein (nonnegative) matrix factorization of histogram matrices. As with most dictionary
learning problems, our objective is not convex but biconvex in the dictionary D and weights Λ and we use ablock-coordinate descent approach. We show that eachof these subproblems can be reduced to an optimization problem involving the Legendre-Fenchel conjugateof the objective, building upon recent work in Cuturi
and Peyré [2016] that shows that the Legendre-Fenchel
conjugate of the entropy regularized Wasserstein distance and its gradient can be obtained in closed form.
We show in Section 4 that these fast algorithms are order of magnitudes faster than those proposed in Sandler and Lindenbaum [2009], whose experiments wereplicate. Finally, we show that the features used
to describe dictionary elements can be different fromthose present in the original histograms. We showcase this property to carry out cross-language semantic
analysis: we learn topics in French using databases ofEnglish texts. A Matlab implementation of our methods and scripts to reproduce the experiment in the introduction are available at http://arolet.github.
io/wasserstein-dictionary-learning/.
Notations If X is a matrix, Xi denotes its ith line,
xj its jth column and Xij its element at the ith line
and jth column. For x, y ∈ Rn , hx, y i is the usual
dot product between x and y. For X, Y ∈ Rn×m ,Pm
def.
hX, Y i = tr(X T Y ) = i=1 hXi , Yi i is the Frobenius
dot-product between matrices X and Y . If A, B areA) denotes
two matrices of the same size, A B (resp. B
the coordinate-wise product (resp. quotient) between
A and B. Σn is the set of n-dimensional histograms:
def. 
Σn = q ∈ Rn+ | hq, 1 i = 1 . If A is a matrix, A+
is its Moore-Penrose pseudoinverse. Exponentials andlogarithms are applied element-wise to matrices andvectors.
If f : Ω ⊂ Rn → R is convex, its Legendre conjugate

Antoine Rolet, Marco Cuturi, Gabriel Peyréf ∗ : x ∈ Rn 7→ maxhx, y i − f (y) is convex. If g : Ω ⊂
y∈Ω

∗

n

R → R is concave, then g∗ = −(−g) is concave andis defined as ∀x ∈ Rn , g∗ (x) = min −hx, y i − g(y).
y∈Ω2Regularized Wasserstein Distances

We start this section by defining Wasserstein distances
for histograms and then introduce their entropic regularization.
2.1

Definition of Wasserstein Distances

Let p ∈ Σn , q ∈ Σs . The polytope of transportationplans between p and q is defined as follows:)(
T1 = pn×s.
U (p, q) = T ∈ R+ s.t.
TT1 = q
Let M ∈ Rn×s
be matrix of costs. The optimal trans+
port cost between p and q with respect to M isdef.

W (p, q) =min hM, T i.

T ∈U (p,q)(1)

two vectors u, v is known as a matrix balancing problem, and is typically solved using Sinkhorn’s [1967]
algorithm. This algorithm has linear convergence, andrequires O(ns) operations at each iteration. The benefits of using an entropy regularization are not justcomputational. In many cases, the linear program defined in Equation (1) does not have a unique solution.
Because of this, W is not differentiable with respect
to either its first or second variable. Wγ , on the other
hand, is differentiable as soon as γ > 0. Problems inwhich we want to optimize an objective depending on
the Wasserstein distance are thus harder to solve with
γ = 0, as argued in Cuturi and Peyré [2016] when trying to compute, for instance, Wasserstein barycenters.
2.3

Legendre Transform

We show in Section 3 that the optimization problems
involved for dictionary learning with a Wasserstein error term can be solved using dual problems whose objectives involve the Legendre-Fenchel conjugate of thesmoothed Wasserstein distance. To abbreviate formulas, we use the following notation for a given p ∈ Σn :def.

Hp = q 7→ Wγ (p, q).
The optimization problem described above is a linear program. It can be seen as a minimum cost network flow problem, and solved with n3 log n complexity when n = s [Orlin, 1993]. In most applications,
M is a pairwise distance matrix. Namely, there exists a metric space (Ω, d) and elements y1 , · · · , yn andz1 , · · · , zs in Ω such that Mij = d(yi , zj ). In that case
W is a distance [Villani, 2009]. We call M the ground
metric of W .

Cuturi and Peyré [2016] showed that the Legendre
transform of the entropy regularized Wasserstein distance, as well as its gradient, can be computed inclosed form:
Hp∗ (g) = γ (E(p) + hp, log Kα i) ,p .
∇Hp∗ (g) = α
KT
Kα
def.

2.2

Entropy Regularization

Solving problems with Wasserstein distance fitting errors can require solving several costly optimal transport problems. As a minimum of affine functions, theWasserstein distance itself is not a smooth functionof its arguments [Cuturi and Doucet, 2014]. To avoid
both of these issues, Cuturi [2013] proposed to smooththe optimal transport problem with an entropic term:def.

Wγ (p, q) =min hM, T i − γh(T ),T ∈U (p,q)(2)

where h is the (strictly concave) entropy function:def.

h(T ) = −hT, log T i.

(3)

3

Wasserstein Dictionary Learning

3.1

Problem Formulation

Let X ∈ (Σn )m be a matrix of m vectors in the nsimplex. Let k be a number of dictionary elements
fixed in advance. We consider the problemminΛ∈Rk×m ,D∈Rs×k

The optimal transport T ? solution of Equation (2) is
unique and can be found by computing two vectorsu ∈ Rn+ , v ∈ Rs+ such that diag(u)K diag(v) ∈ U (p, q),
where K = e−M/γ . The optimal solution is then thenT ∗ = diag(u)K diag(v). The problem of finding thesedef.

where K = e−M/γ and α = eg/γ . Moreover, they
showed that for p ∈ Σn and g ∈ Rs , ∇Hp∗ (g) ∈ Σsand that ∇Hp∗ is γ1 -Lipschitz. When γ = 0, Hp∗ is notdifferentiable anymore but elements of its subgradient can be computed efficiently, notably for Euclidean
point clouds [Carlier et al., 2014].

m
X

Hxi (Dλi ) + R(Λ, D)

(4)

i=1

s.t. DΛ ∈ Σms.

Problem (4) is convex separately (but not jointly) in
D and Λ as long as R is convex. We propose in what
follows to use a block-coordinate descent on D and Λ.

Fast Dictionary Learning with a Smoothed Wasserstein Loss

Sandler and Lindenbaum [2009] show that when R =0, γ = 0 and either D or Λ is fixed, Equation (4)is a linear program of dimensions (m × n × s) withm×(n×s+n+s) constraints, each involving 1, n or t×k
variables). Representing these constraints is challenging for common sized datasets, and solving such problems is usually intractable. They proposed to replace
the Wasserstein distance by an approximation [Shirdhonkar and Jacobs, 2008], for which the gradients areeasier to compute. However this approximation can
only be used when M is a distance matrix in a Euclidean space of small dimension. We propose to consider instead a positive regularization strength γ > 0.
This allows us to consider any cost matrix M , rather
than only pairwise distance matrices, and makes the
optimization problems smooth and better behaved, inthe sense that when D or Λ are full rank, the optimizers of each block update is unique. We propose nextin §3.3 an entropy regularization on the columns of Dand Λ to enforce positivity of these coefficients.
3.2

Wasserstein Dictionary Learning

Mixture Weights Update We consider here thecase where the dictionary D is fixed, and our goal is
to compute mixture weights Λ

argminmX

Hxi (Dλi ),Λ∈Rk×m i=1

s.t. DΛ ∈ Σm
s .

(5)

This problem can be solved using a gradient descent,but computing the gradient is equivalent to evaluating each Hxi (Dλi ) for i = 1, · · · , m, that is solving mintermediate matrix scaling problems. We propose touse duality instead and attack the problem by exploiting the fact that Hx∗i have closed form gradients.
Theorem 1. Let Λ? be a solution of Problem (5). Λ?
satisfies Dλ?i = ∇Hx∗i (gi? ) for i = 1, . . . , m, withgi? ∈ argminHx∗i (g) s.t. DT g = 0.

straints, so strong duality holds and its dual reads
maxG∈Rs×m Λ∈Rk×m ,Q∈Σms

G∈Rs×m Λ∈Rk×m

Λ∈Rk×mG∈R

Λ∈Rk×m

G∈R

m
X

m
X

Hxi (qi ) − hqi , gi ii=1

Hx∗i (gi )

(7)

i=1
m
X

= max
min
s×m

If DT G 6= 0, thenhΛi , DT gi i − Hx∗i (gi ).

i=1

min

m
P

hλi , DT gi i = −∞. SinceΛ∈Rk×m i=1

Hx∗i (gi ) is finite for all i, the maximum over G is realized only if DT G = 0. The problem becomesmax

G∈Rs×ms.t. DT G=0−m
X

Hx∗i (gi ) =

i=1

m
X

maxs

g∈R
i=1
s.t. DT g=0−Hx∗i (g).

This is the same optimization problem as in Equation (6) and it has only one solution G? . Them first
order conditions of (7) are DΛ? = ∇Hx∗i (gi? ) i=1 . If
D is full rank, this linear equation has a unique solution.
Remark 1. Here, for i = 1 . . . m, Dλ?i is in the simplex because ∇Hx∗i (gi ) is itself in the simplex (see section 2.3). However the column of Λ? are not required
to be in the simplex and could even possibly take negative values. If all the columns of D are in the simplex,
then, however, columns of Λ? need to sum to 1.
We solve Equation (6) with a projected gradient descent and then recover Λ? by solving the linear equamtion DΛ? = ∇Hx∗i (gi? ) i=1 .
Dictionary Update Assuming weights Λ are fixed,
our goal is now to learn the dictionary matrix D.
Theorem 2. Let D? be a solution ofD∈Rs×kProof. Let us introduce the variable Q = DΛ. Problem (5) becomes

Q∈Σt

= max
min hDΛ, G i −
s×m

min

Moreover if D is full-rank this solution is unique.

Hxi (qi ) + hDΛ − Q, G i

i=1

min hDΛ, G i + minm

= max

(6)

g∈Rs

m
X

min

m
X

Hxi (Dλi ) s.t. DΛ ∈ Σm
s .

i=1

m
D? satisfies D? Λ = ∇Hx∗i (gi? ) i=1 , with
G? ∈ argminm
X

G∈Rs×m i=1

Hx∗i (gi ) s.t. GΛT = 0.

(8)Moreover if Λ is full-rank this solution is unique.
minΛ∈Rk×m ,Q∈Σmsm
X

Hxi (qi ) s.t. DΛ = Q.

i=1

It is a convex optimization problem with affine con-

The proof is similar to that of Theorem 1. We solve
Equation (8) with a projected gradient descent and?
?
then recover D by solving the linear equation D Λ =∗? m∇Hxi (gi ) i=1 .

Antoine Rolet, Marco Cuturi, Gabriel Peyré

3.3

Wasserstein NMF

In order to enforce non-negativity constraints on thevariables, we consider the problem
min

m
X

Λ∈Σm
k
i=1
D∈Σk
s

Hxi (Dλi ) − ρ1 E(Λ) − ρ2 E(D)3.4,

where E is defined for matrices which columns are in
the simplex as E(A) = hA, log A i. This regularization allows us to derive similar results as those in Section 3.2 which we can use to find non-negative iterates
for D and Λ efficiently.
Enforcing Positive Weights Problem (9) with afixed dictionary is convex but as in Section 3.2 thegradient of the objective is computationally expensive.
We have a similar duality result:
Theorem 3. The solution Λ? of
min

Λ∈Σm
k

Hxi (Dλi ) − ρ1 E(λi ) s.t. DΛ ∈ Σm
s ,

i=1
e−D

is such that ∀i = 1 . . . m, λ?i =T g ? /ρ1i

T ?
he−D gi /ρ1 ,1 i

T
k
X
m
e−GΛi /ρ2 Λi∗.
∇Hxi (gi ) i=1 −T
he−GΛi /ρ2 , 1 i
i=1

(9)

s.t. DΛ ∈ Σmsm
X

The proof is similar to that of Theorem 3. We solve
Equation (11) with an accelerated gradient scheme.
The gradient of the objective of in Equation (11) is, withgi? ∈ argmin Hx∗i (g) − ρ1 E∗ (−DT g/ρ1 ).

(10)

Convergence

As pointed by Sandler and Lindenbaum [2009], thealternate optimization process generates a sequenceof lower bounded non-increasing values for the objective of Problem (4), so the sequence of objectives converges. When, moreover, we use an entropy regularization (ρ1 , ρ2 > 0, §3.3), successive updates for D andΛ remain in the simplex, which is compact, and thus
satisfy the conditions of [Tropp, 2003, Theorem 3.1],
taking into account that the hypothesis made in thattheorem that the divergence is definite is not actuallyused in the proof. Thus every accumulation point of
the sequences of iterates of D and Λ is a generalized
fixed point. Moreover, if the iterates remain of fullrank, then Theorem 3.2 in the same reference applies,
and the sequences either converge or have a continuum of accumulation points. Although this full rank
hypothesis is not guaranteed to hold, we observe that
it holds in practice when the entropy regularization
term does not dominate the objective.

g∈Rs

The proof, similar to that of Theorem 1, is given in appendix. The objective and gradient of the optimization problem in Equation (10) can be computed inclosed form with the following formulas:E∗ (x) = − loghex , 1 i,

∇E∗ (x) = −ex.
hex , 1 i

We solve Equation (10) with an accelerated gradientscheme [Nesterov, 1983]. The gradient of the objective
in Equation (10) is!mT
e−D gi /ρ1∗∇Hxi (gi ) − D −DT g /ρ.
i
1, 1 ihe
i=1

Enforcing Positive Dictionaries
have the following theorem:
Theorem 4. The solution D? of
min

D∈Σk
s

is d?i =m
X

Hxi (Dλi ) − ρ2

i=1

E(di ) s.t. DΛ ∈ Σm
s ,

i=1

? T
e−G Λi /ρ2−G? ΛT /ρ2
i
he
,1 i

G? ∈ argminkX

Similarly, wemXG∈Rs×m i=1, withHx∗i (gi ) −

k
X
i=1

ρE∗ (−GΛTi /ρ2 ). (11)3.5

Implementation

Projection Step for Unconstrained Dictionary
Learning We solve Equations (6) and (8) with projected gradient descent methods. The orthogonal projector of the optimization problem is projKer(DT ) :=G 7→ G − DD+ G in Equation (6) and projKer(Λ) :=G 7→ G − GΛ+ Λ in Equation (8). Precomputing DD+(resp. Λ+ Λ) uses O(s2 ) (resp. O(m2 )) memory space,
and then the projection is performed in complexity
O(s2 × m) (resp. O(s × m2 )). When either s or m islarge, storing such a matrix is too expensive and leads
to slowdowns due to memory management. In such acase, we can precompute D+ (resp. Λ+ ), which takesO(s×k) (resp. O(m×k)) memory space, and computeprojKer(DT ) (G) as G − D(D+ G) (resp. projKer(Λ) (G)as G − (GΛ+ )Λ) in O(s2 × m2 × k 2 ) operations.
Parallelization of the Dictionary Update Parallelization on multiple processes is easy for the weights
updates because each weight vector λi can be computed independently. The dictionary updates however cannot be reduced to completely independent subproblems. Indeed the constraint in Equation (8) makes
a dependence on the columns of D. Similarly the gradient of the objective in Equation (11) cannot be separated into independent sub-problems.

Fast Dictionary Learning with a Smoothed Wasserstein Loss

We show how to use parallel processes to speed-up theunconstrained dictionary updates. The most computationally expensive part it to solve the optimization
problem of Equation (8). The objective and gradient of this problem can be computed independently
for each column. Then we can gather the gradient ona single process and project it. Since the constraint
is linear we can project directly the gradient before
computing the step-size of the descent, so that if thiscomputation involves computing the objective (like abacktracking line-search does for example) the projection does not need to be repeated.
We also propose a scheme to partially parallelize thepositive dictionary updates. The objective and gradient of the optimization problem in Equation (11)T
are found by computing e−GΛ , which cannot be computed separately on columns of G. An efficient way toTstill compute e−GΛ in parallel is to split G columnwise into (G(1) , . . . , G(p) ) where p is the number of(i) T
processes available, and compute e−G Λ on proT
cess i. The managing process computes e−GΛ as
Qp
(i) T
−G Λ
(here the product is point-wise) andi=1 egives the result to all the other processes so thatthey can finish computing the gradient. By doing so
most of the work is done in parallel and each process
only shares a matrix of size s × k twice per gradient/objective calculation. Since usually k  m thisallows to use all the available processes while keeping
communication overhead low.

4

Experiments

4.1

Face Recognition

We reproduce here the face recognition experimentof Sandler and Lindenbaum [2009] on the ORL
dataset [Samaria and Harter, 1994] with the same
preprocessing, classification and evaluation method inorder to compare computation time. Each image isdownsampled so that its longer side is 32. We represent images as column vectors that we normalize
so that they sum to 1 and store them in matrix X.
The cost matrix M is the Euclidean distance between
pixels. For evaluation, the dataset is split evenly in
two, trained on one set and tested on the other several
times, and we take the take best classification performance obtained. Table 1 shows the classification accuracy obtained with unconstrained Wasserstein Dictionary Learning (Section 3.2). The results are comparable to those of Sandler and Lindenbaum [2009].
Learning the dictionary and coefficients with a Matlab
implementation of our algorithm on an single core of
a 2.4Ghz Intel Quad core i7 CPU with k = 40 takes

Table 1: Classification accuracy for the face recognition task on the ORL dataset.
k
10
20
30
40
50
γ = 1/30
93% 95.5% 97% 96.5% 96%
γ = 1/50
91%
95% 95% 97% 94.5%
Sandler09 94.5% 90.5% 95% 96.5% 97%

on average 20s for γ = 1/30 and 90s for γ = 1/50,
while Sandler and Lindenbaum [2009] report up to 20
minutes just for the D step with a comparable CPU.
The whole NMF can take up to 10 minutes when we
use the entropy positivity barrier with ρ1 = ρ2 = 1/10.

4.2

Semantic Analysis

The goal of semantic analysis is to extract a few representative histograms of words (a.k.a. topics) fromlarge corpora of texts. NMF is a classical approach tocarry out this task, either as is [Lee and Seung, 1999],
or with a different factorization scheme to model thedata generation process, such as Probabilistic Latent
Semantic Analysis [Hofmann, 1999], which is a nonnegative factorization of the form X = DΣΛ.
We use the same approach as Lee and Seung [1999] tolearn topics from a database of text with NMF. The
input data is a bag-of-words representation of the documents. Let Y = (y1 , . . . , yn ) be the vocabulary ofthe database, a text document is represented as vectorof word frequencies: Xij is the frequency of the word
yi in the j th text. We then learn a factorization DΛwith NMF. The topics that we learn are the columnsof D. The cost of the factorization is usually its Euclidean distance or Kullback-Leibler divergence to X.
In order to use a Wasserstein cost instead, we needa meaningful cost for transporting words from one toanother.
Recent works [Pennington et al., 2014, Zou et al.,2013], building upon earlier references [Bengio et al.,2003], propose to compute Euclidean embeddings forwords such that the Euclidean or cosine distances between the respective image of two words corresponds
to some form of semantic discrepancy between thesewords. As recently shown by Kusner et al. [2015],
these embeddings can be used to compare texts using
the toolbox of optimal transport: Bag-of-words histograms can be compared with Wasserstein distances
using the Euclidean metric between the words as theground metric M . We leverage these results to learn
topics from a text database using Wasserstein NMF.

Antoine Rolet, Marco Cuturi, Gabriel Peyré

4.2.1

Datasets

We learned topics on two labeled datasets. For eachdataset, m is the number of documents, n is the vocabulary size and V is the number of labels. (i) BBCsport [Greene and Cunningham, 2006] is a dataset ofnews articles about sports, labeled according to whichsport the article is about, in which we removed stopwords (n = 12, 669, m = 737, V = 5). We split thedataset as a 80/20 training / testing set for classification. (ii) Reuters is a dataset of news articles
labelled according to their area of interest. We usedthe version described in Cardoso-Cachopo [2007], withthe same train-test split for classification, and removed
stop-words and words that appeared only once acrossthe whole corpus (n = 13, 038, m = 7, 674, V = 8).
4.2.2

Monolingual Semantic Analysis

We used a pretrained Glove word embedding [Pennington et al., 2014] to map words to a Euclidean space ofdimension 300. Let g denote this embedding, and Z =(z1 , . . . , zs ) the target vocabulary, that is the words
that are allowed to appear in the topics. We define thecost matrix of the Wasserstein distance as the cosinehg(y ),g(z ) idistance in the embedding: Mij = 1 − kg(yii)kkg(zj j )k .
We then find D and Λ with Wasserstein NMF (WNMF, Section 3.3).
Figure 2 shows a word cloud representation (wordle.
net) of 4 relevant topics for the dataset BBCsport. Depending on the parameters, the full Wasserstein NMF
computation takes from 20 minutes to an hour for
BBCsport and around 10 hours for Reuters using a
Matlab implementation running on a single GPU ofan Nvidia Tesla K80 card.
Target Words Selection Since we can choose as
target words any word that is defined for the embedding, we need a way to select which to use. We
chose to use a list of 3, 000 frequent words in English1 .
Other approaches can be considered such as using thedataset’s vocabulary, tokenized or not, or taking themost frequent words for each class in the dataset.
4.2.3

Figure 2: Word clouds representing 4 of the 15 topics learned on BBCsport in English. Top-left topic:
competitions. Top-right: time. Bottom-left: soccer
actions. Bottom-right: drugs.

Cross-language Semantic Analysis

Lauly et al. [2014] propose a bilingual word representation that maps words in two different languages to thesame Euclidean space. By setting Z, the vocabularyof the topics, as a subset of the words in the target language, we can learn topics in this language. Figure 3
illustrates what we would expect with k = 1, which is
the Wasserstein iso-barycenter problem. We use a pre1
Available at https://simple.wiktionary.org/wiki/
Wiktionary:BNC_spoken_freq

trained embeddings of dimension 40 from Lauly et al.
[2014] in order to learn topics in French. Note that
this method could also learn topics in one language
from a bilingual dataset, or in both languages.
As in Section 4.2.2, we use the cosine distance inthe embedding as the ground metric. Table 4 shows
word cloud representations of 4 relevant topics for thedataset Reuters. Computation times are similar to
those with a target vocabulary in English.

Figure 3: The Wasserstein iso-barycenter of two English sentences with a target vocabulary in French. Arrows represent the optimal transport plan from a text
to the barycenter. The barycenter is supported on
the bold red words which are pointed by arrows. The
barycenter is not equidistant to the extreme points because the set of possible features is discrete.

Target words selection We chose as the target dictionary a list of of 6, 000 frequent words in French2 .
2
Available at http://wortschatz.uni-leipzig.de/Papers/top10000fr.txt

Fast Dictionary Learning with a Smoothed Wasserstein Loss

Table 3: Confusion matrices for BBCsports for k-NN
with W-NMF. Columns represent the ground truth
and lines predicted labels. Labels: athletism (a),cricket (c), football (f), rugby (r) and tennis (t).
English target vocabulary

French target vocabularya c f r t
a 21 0 0 0 0
c 0 25 0 0 0
f 0 0 50 4 1
r 0 0 3 26 0
t 0 0 0 0 19

a c f r t
a 18 0 1 0 2
c 0 22 3 0 2
f 3 0 44 5 6
r 0 3 3 25 1
t 0 0 2 0 9

Figure 4: Word clouds representing 4 of the 24 topics
learned on Reuters in French. Top-left topic: international trade. Top-right: oil and other resources.
Bottom-left: banking. Bottom-right: management
and funding.
Table 2: Text classification error. W-NMFf is the classifier using W-NMF with a French target vocabulary.

Table 4: 10 French words closest to some English
words according to the ground metricfootball
MethodKL-NMF E-NMF W-NMF W-NMFf
Reuters
6.9%
8.2%
6.0%
9.8%
BBCsport
9.4%
12.8%
5.4%
20.8%

bank
tennis4.2.4

Classification Performance

We compared the classification error obtained on thetwo datasets with our method to those obtained byusing the coefficients produced by Euclidean NMF (ENMF) and Kullback-Liebler NMF (KL-NMF). The
classifier we used is k-NN using the Hellinger distance,with k selected by 10-fold cross-validation on the training set, using the same partitions for all methods. We
set the number of topics to 3V. Parameters γ, ρ1 andρ2 were set to be as small as we could (small values
can make the gradients infinite because of machine precision) without a particular selection procedure. See
supplementary materials for a representation of all the
topics of every method.
Wasserstein NMF with a target vocabulary in English
performs better on this auxiliary task than Euclidean
or KL NMF. Although this does not prove that the
topics are of better quality, it shows that Wasserstein
NMF can drastically reduce the vocabulary size without loosing discrimitative power. As we can see in
Figures 2, 4, the topics themselves are semanticallycoherent and related to the datasets’ content.
The classification error for W-NMF with a French target vocabulary on BBCsports is rather bad, although
the topics are coherent and related to the content of
the articles. The confusion matrix (Table 3) showsthat more than half of the articles about tennis are
misclassified. In fact, the other methods produce a
topic about tennis, but W-NMF with a French dictionary does not. Table 4 shows the French words

football supporters championnat sportives
sportifs joueurs sportif jeux matches sportbanque banques bei bancaire federal bankemprunts reserve crédit bancaires
bienfaiteurs murray ex-membre balletb92 sally sylvia markovic hakim socialocommuniste

closest to some English query words according to theground metric. While the closest words to football andbank are semantically related to their query word, theclosest words to tennis are not. This illustrates how
our method relies on the ground metric, given by wordembeddings in this case.

5

Conclusion

We show how to efficiently perform dictionary leaning and NMF using optimal transport as the data fitting term, with an optional entropy positivity barrier.
Our method can be applied to large datasets in highdimensions and does not require any assumption on
the cost matrix. We also show that with this data
fitting term, the reconstruction DΛ can use different
features than the data X. Other than our application
to cross-language semantic analysis, this can be used
for example to reduce the number of target featuresby quantization for the dictionary while keeping theoriginal features for the dataset.
While we only consider entropy as a barrier for positivity in this work, our approach can be generalized
to other regularizers, as long as the gradient of R? or
its proximal operator can be computed efficiently. We
believe that extensions to other classes of regularizers
is an interesting area for future work.

Antoine Rolet, Marco Cuturi, Gabriel Peyré

References
Y. Bengio, R. Ducharme, P. Vincent, and C. Janvin.
A neural probabilistic language model. The Journalof Machine Learning Research, 3:1137–1155, 2003.
A. Cardoso-Cachopo. Improving Methods for Singlelabel Text Categorization. PdD Thesis, Instituto
Superior Tecnico, Universidade Tecnica de Lisboa,
2007.
G. Carlier, A. Oberman, and E. Oudet. Numericalmethods for matching for teams and wassersteinbarycenters.
arXiv preprint arXiv:1411.3602v1,
2014.
A. Chambolle and C. Dossal.
On the convergence of the iterates of the fast iterative shrinkage/thresholding algorithm. Journal of Optimization Theory and Applications, 166(3):968–982, 2015.
M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems 26, pages 2292–2300,
2013.
M. Cuturi and D. Avis. Ground metric learning. The
Journal of Machine Learning Research, 15(1):533–
564, 2014.
M. Cuturi and A. Doucet. Fast computation of wasserstein barycenters. In 31st International Conference
on Machine Learning (ICML-14), 2014.

D.D. Lee and H.S. Seung. Algorithms for non-negative
matrix factorization. In Advances in neural information processing systems, pages 556–562, 2001.
Y. Nesterov. A method of solving a convex programming problem with convergence rate o (1/k2). In
Soviet Mathematics Doklady. Vol. 27. No. 2., 1983.
James B Orlin. A faster strongly polynomial minimum
cost flow algorithm. Operations research, 41(2):338–
350, 1993.
P. Paatero and U. Tapper. Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values. Environmetrics, 5(2):111–126, 1994.
J. Pennington, R. Socher, and C.D. Manning. Glove:
Global vectors for word representation. Proceedingsof the Empiricial Methods in Natural Language Processing (EMNLP 2014), 12:1532–1543, 2014.
J. Rabin and N. Papadakis. Convex color image segmentation with optimal transport distances. In
Scale Space and Variational Methods in Computer
Vision, pages 256–269. Springer, 2015.
Y. Rubner, C. Tomasi, and L.J. Guibas. A metric for
distributions with applications to image databases.
In Computer Vision, 1998. Sixth International Conference on, pages 59–66. IEEE, 1998.

M. Cuturi and G. Peyré. A smoothed dual approachfor variational Wasserstein problems. SIAM Journal
on Imaging Sciences, 2016. to appear.

F.S. Samaria and A.C. Harter. Parameterisation of astochastic model for human face identification. In
Applications of Computer Vision, 1994., Proceedings of the Second IEEE Workshop on, pages 138–
142. IEEE, 1994.

D. Greene and P. Cunningham. Practical solutions to
the problem of diagonal dominance in kernel document clustering. In Proc. 23rd International Conference on Machine learning (ICML’06), pages 377–
384. ACM Press, 2006.

R. Sandler and M. Lindenbaum. Nonnegative matrix
factorization with earth mover’s distance metric. In
Computer Vision and Pattern Recognition, 2009.
CVPR 2009. IEEE Conference on, pages 1873–1880.
IEEE, 2009.

T. Hofmann. Probabilistic latent semantic analysis. In
Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, UAI’99, pages 289–
296, San Francisco, CA, USA, 1999. Morgan Kaufmann Publishers Inc. ISBN 1-55860-614-9.

M.N. Schmidt, J. Larsen, and F.T. Hsiao. Wind noise
reduction using non-negative sparse coding. In Machine Learning for Signal Processing, 2007 IEEE
Workshop on, pages 431–436. IEEE, 2007.

M.J. Kusner, Y. Sun, N.I. Kolkin, and K.U. Weinberger. From word embeddings to document distances. In Proceedings of The 32nd International
Conference on Machine Learning, 2015.
S. Lauly, H. Larochelle, M. Khapra, B. Ravindran,
V.C. Raykar, and A. Saha. An autoencoder approach to learning bilingual word representations.
In Advances in Neural Information Processing Systems, pages 1853–1861, 2014.
D.D. Lee and H.S. Seung. Learning the parts of objects
by non-negative matrix factorization. Nature, 401(6755):788–791, 1999.

S. Shirdhonkar and D.W. Jacobs. Approximate earth
movers distance in linear time. In Computer Vision
and Pattern Recognition, 2008. CVPR 2008. IEEE
Conference on, pages 1–8. IEEE, 2008.
R. Sinkhorn. Diagonal equivalence to matrices withprescribed row and column sums. The American
Mathematical Monthly, 74(4):402–405, 1967.
J. Solomon, F. de Goes, Pixar Animation Studios,
G. Peyré, M. Cuturi, A. Butscher, A. Nguyen,
T. Du, and L.J. Guibas. Convolutional wasserstein
distances: Efficient optimal transportation on geometric domains. ACM Transactions on Graphics(Proc. SIGGRAPH 2015), 2015.

Fast Dictionary Learning with a Smoothed Wasserstein Loss

J.A. Tropp. An alternating minimization algorithm for
non-negative matrix approximation, 2003.
C. Villani. Optimal transport: old and new, volume338. Springer Verlag, 2009.
G. Zen, E. Ricci, and N. Sebe. Simultaneous ground
metric learning and matrix factorization with earth
mover’s distance. In Pattern Recognition (ICPR),
2014 22nd International Conference on, pages 3690–
3695. IEEE, 2014.
S. Zhang, W. Wang, J. Ford, and F. Makedon. Learning from incomplete ratings using non-negative matrix factorization. In SDM, volume 6, pages 548–552.
SIAM, 2006.
W.Y. Zou, R. Socher, D.M. Cer, and C.D. Manning.
Bilingual word embeddings for phrase-based machine translation. In EMNLP, pages 1393–1398,2013.