Warpspeed Computation of Optimal Transport, Graph Distances, and Embedding Alignment
Authors: Johannes Klicpera
Marten Lienen
Stephan Günnemann
Venue:  (2021)
Keywords: Optimal transport
sinkhorn distance
locality sensitive hashing
nyström method
graph neural networks
embedding alignment

Abstract: Optimal transport (OT) is a cornerstone of many machine learning tasks. The current best practice for computing OT is via entropy regularization and Sinkhorn iterations. This algorithm runs in quadratic time and requires calculating the full pairwise cost matrix, which is prohibitively expensive for large sets of objects. To alleviate this limitation we propose to instead use a sparse approximation of the cost matrix based on locality sensitive hashing (LSH). Moreover, we fuse this sparse approximation with the Nyström method, resulting in the locally corrected Nyström method (LCN). These approximations enable general log-linear time algorithms for entropy-regularized OT that perform well even in complex, high-dimensional spaces. We thoroughly demonstrate these advantages via a theoretical analysis and by evaluating multiple approximations both directly and as a component of two real-world models. Using approximate Sinkhorn for unsupervised word embedding alignment enables us to train the model full-batch in a fraction of the time while improving upon the original on average by 3.1 percentage points without any model changes. For graph distance regression we propose the graph transport network (GTN), which combines graph neural networks (GNNs) with enhanced Sinkhorn and outcompetes previous models by 48%. LCN-Sinkhorn enables GTN to achieve this while still scaling log-linearly in the number of nodes.

Under review as a conference paper at ICLR 2021
WARPSPEED COMPUTATION OF OPTIMAL TRANSPORT, GRAPH DISTANCES, AND EMBEDDING ALIGNMENT
Anonymous authors Paper under double-blind review
ABSTRACT
Optimal transport (OT) is a cornerstone of many machine learning tasks. The current best practice for computing OT is via entropy regularization and Sinkhorn iterations. This algorithm runs in quadratic time and requires calculating the full pairwise cost matrix, which is prohibitively expensive for large sets of objects. To alleviate this limitation we propose to instead use a sparse approximation of the cost matrix based on locality sensitive hashing (LSH). Moreover, we fuse this sparse approximation with the Nyström method, resulting the locally corrected Nyström method (LCN). These approximations enable general log- -ig-linear time algorithms for catropy-regllariee OT that perform well even in complex. high- sigh--immnsional spaces. We thoroughly demonstrate these advantages via theoretical analysis and by evaluating multiple approximations both directly and as a component of two real-world models. Using approximate Sinkhorn for unsupervised word embedding alignment enables us to train the model full-batch in a fraction of the time while improving upon the original on average by 3.1 percentage points without any model changes. For graph distance regression we propose the graph transport network (GTN), which combines graph neural networks (GNNs) with enhanced Sinkhorn and outcompetes previous models by 48 %. LCN-Sinkhorn enables GTN to achieve this while still scaling log-linearly in the number of nodes.
INTRODUCTION
Measuring the distance between two distributions or sets of objects is central problem in machine learning. One common method of solving this is optimal transport (OT). OT is concerned with the problem of finding the transport plan for moving source distribution (e.g. pile of earth) to sink distribution (e.g. a construction pit) with the cheapest cost w.r.t. some pointwise cost function (e.g. the Euclidean distance). The advantages of this method have been shown numerous times, e.g. in generative modelling (Arjovsky et al. 2017 Bousquet et al. 2017 Genevay et al., 2018), loss functions (Frogner et al. 2013)) set matching (Wang al. 2019). or domain adaptation (Courty et al.) 2017). Motivated by this, many different methods for accelerating OT have been proposed in recent years (Indyk & Thaper 2003) Paapaakii al. 20144] Backurs et al., 2020). However, most of these approaches are specialized methods that do not generalize to modern deep learning models, which rely on dynamically changing nigh-dimensional embeddings.
In this work we aim to make OT computation for point sets more scalable by proposing two fast and accurate approximations of entropy-regularize optimal transport: Sparse Sinkhorn and LCN- Sinkhorn, the latter relying on our newly proposed locally corrected Nyström (LCN) method. Sparse Sinkhorn uses a sparse cost matrix to leverage the fact that in entropy- regulapieee OT (also known as the Sinkhorn distance) (Cuturi 2013) often only each point's nearest neighbors influence the result. LCN-Sinkhorn extends this approach by leveraging I.CN, a general similarity matrix approximation that fuses local (sparse) and global (low-rank) approximations, allowing to simultancously capture both kinds of behavior. LCN-Sinkhoo Sinkhorn thus fuses sparse Sinkhorn and Nyström-Sinkhomm (Altschuler et al. 2019). Both sparse Sinkhorn and LCC--innhon run in log- time.
We theoretically analyze these approximations and show that sparse corrections can lead to significant improvements over the Nyström approximation. We furthermore validate these approximations by showing that they are able to reproduce both the Sinkhorn distance and transport plan significantly better than previous methods across wide range of regularization parameters and computational

---

Under review as a conference paper at ICLR 2021
Multiscale OT
Nyström Sinkhorn
Sparse Sinkhorn
LCN-Sinkhorn
(porfet correlatio)
Diagonal
PNys.5
PLCNN
pMs
Pa
P
Pij
Po
Pg
Pij
Figure I: The proposed methods (sparse and LCN-Sinkhorn) show a clear correlation with the full Sinkhorn transport plan, opposed to previous methods. Entries of approximations (y-axis) and full Sinkhorn (X-axis) for pre-aligned word embeddings (EN-DE). Color denotes sample density.
budgets (as e.g. demonstrated in Fig. We then show the impact of these improvements by employing Sinkhorn approximations end-to-end in two high-impact machine learning tasks. First, we incorporate them into Wasserstein Procrustes for word embedding alignment (Grave al. 2019). LCN-Sinkhorn improves upon the original method's accuracy by 3.1 percentage points using a third of the training time without any further model changes. Second, (GTN), which combines graph neural networks (GNNs) with optimal transport, and further improve via learnable unbalanced OT and multi-head OT. GTN with LCN-Sinkhorn is the first model that both overcomes the bottleneck of using a single embedding per graph and scales log-linearly in the number of nodes. In summary, our paper's main contributions are:
develop the graph transport network
Locally Corrected Nyström (LCN), flexible. log-linear time approximation for similarity matrices, leveraging both local (sparse) and global (low-rank) approximations.
Eegroiy-regullrized optimal transport (a.k.a. Sinkhorn distance) with log-linear runtime via sparse Sinkhorn and CN--inkhomn. These are the first tog-lineau approximations that are stable enough to substitute full entropy regullareed OT in models that leverage high nigh-dimensional spaces.
The graph transport network (GTN), which combines graph neural network (GNN) with multi- head unbalanced LCN- CCN-Sinkhorn. GTN both sets the state of the art on graph distance regression and still scales log- oineeilyryy in the number of nodes.
2 SPARSE SINKHORN
Entropy-reguulriized optimal transport. In this work we focus on optimal transport between two discrete sets of points. We furthermore add entropy regularization, which enables fast computation and often performs better than regular OT (Cuturi) 2013). Formally, given two categorical distributions modelled via the vectors p E R" and q E Rnth supported on two sets of points Xp {_pp...... Tpn
and Xg {Teg...... yym} in gd and the cost function c Rd R + E (e.g. the squared La
distance) giving rise to the cost matrix Cty c(Epi. Dgz) we aim find the Sinkhorn distance d and the associated optimal transport plan (Cuturi 2013)
d=
min P. C)r (H)
s.t. Plm p. p' 1, q:
(1)
P
with the Frobenius inner product (., .)F and the entropy
Em Pij log Poj. Note 1/2 that d includes the entropy and can thus be negative, while Cuturi (013)) originally used daturi, (P, C. This optimization problem can be solved by finding the vectors $ and that normalize the columns and rows of the matrix P diag(3 K diag(t) with the similarity matrix K,, = so that P1 m p and pT = q. This is usually achieved via the Sinkhorn algorithm, which initializes Im and then updates them alternatingly via
the normalization vectors (1) 1,, and (1))
@ (Kt(i-)),
(00 2 (K" (()
until convergence, where @ denotes elementwise division.
Sparse Sinkhorn. The Sinkhorn algorithm is faster than non-regularized EMD algorithms, which run in On2m log n log(n max(C))) TTarjan 1997). However, its computational cost is still quadratic
2

---

Under review as a conference paper at ICLR 2021
in time, i.e. Onm)), which is prohibitively expensive for large n and m. We propose overcome this by observing that the matrix K, and hence also P. is negligibly small everywhere except at each point's 's closest neighbors because of the exponential used in K's computation. We propose to leverage this by approximating C' via the sparse matrix C where
Cij 0c
if Tpi and Faj are "near", otherwise.
(3)
KdP and P' follow according to the definitions of K and P. In this work we primarily consider neighbors with distance lower than r1 "near". Finding such neighbors can be efficiently solved via locality sensitive hashing (LSH) on Xp U U.y.
Locality sensitive hashing. LSH tries to filter "near" from "far" data points by putting them into different hash buckets. Points closer than a certain distance "1 are put into the same bucket with probability at least pt, while those beyond some distance r2 c ri with e are put into the same bucket with probability at most & Pl. There is a plethora of LSH methods for different cost functions (Wang et al., 2014 Shrivastava Li, 2014). so we do not have to restrict our approach to a limited set of functions. In this work we focus on cross-polytope LSH (Andoni et al.) 2015) and k-means LSH (Paulevé et al. 2010), depending on the cost function (see App. H). Sparse Sinkhorn with LSH scales log-linearly with the number of points, i.e. On log n) for R m (see and App. K for details). Unfortunately, LSH can fail when e.g. the cost between pairs very similar (see App. B) However, we can alleviate these limitations by fusing K'P with the Nyström approximation.
3 LOCALLY CORRECTED NYSTRÖM AND
Nyström method. The Nyström method is popular way of approximating similarity matrices that provides performance guarantees for many important tasks (Williams & Seeger [2001] Musco & Musco (017). It approximates a positive semi-definite (PSD) similarity matrix K via its low-rank decomposition KNys UA -V.. Since the optimal decomposition via SVD is too expensive to compute, Nyström instead chooses a set of landmarks I, {.11.-.. my and obtains the matrices
via U., k( (ty). A.j 1 k( Diy)) and V12 = k( Fgj), where k(x11,2) 222 is an arbitrary
for Sinkhorn. Common methods of choosing landmarks from X, U Xg are uniform and ridge leverage score (RI-) sampling. We instead focus on k-means Nyström and sampling via k-means++, which we found be significantly faster than recursive RLS sampling Zhhang et al. 2008) and perform better than both uniform and RLS sampling (see App. H).
PSID kernel, e.g. k(2r1, 272) = e stelmes
Sparse vs. Nyström. Exponential kernels like the one used for K (e.g. the Gaussian kernel) typically have a reproducing kernel Hilbert space that is infinitely dimensional. The resulting Gram matrix K thus always has full rank. A low-rank approximation like the Nyström method can therefore only account for its global structure and not the local structure around each point æ. As such, it is ill ill-suited for any moderately low entropy regularization parameter, where the transport matrix P resembles permutation matrix. Sparse Sinkhorn, on the other hand, cannot account for global structure and instead approximates all non-selected distances as infinity. It will hence fail if more than a handful of neighbors are required per point. These approximations are thus opposites of each other, and such not competing but rather complementary approaches.
Locally corrected Nyström. Since we know that the entries in our sparse approximation exact, fusing this matrix with the Nyström method is rather straightforward. For all non non-zero zero values in the sparse approximation K we first calculate the corresponding Nyström approximations, obtaining the sparse matrix K To obtain the locally corrected Nyström (LCN) approximation we remove Nys these entries from K Kyys and replace them with their exact values, i.c.
KICN = KNys + K = KNYs K K'P
(4)
LCN-Sinkhorn. To obtain the approximate transport plan PL.N we run the Sinkhorn algorithm with KLCN instead of K. However, we never fully instantiate KLCN. Instead, we only save the decomposition and directly use these parts in Eq. (2) via Kiext UA ---V) + K similarly
to Altschuler et al. (2019). As result we obtain the decomposition of PLAN Pkys + P Po Pw + psp_ P and the approximate distance (using Lemma A from Altschuler et al. 2019
d deN.e A Pu Pw 1m + Pr + 8 111, + 1 PTt).
(5)
3

---

Under review as a conference paper at ICLR 2021
This approximation scales log- log-linearly with dataset size (see Ssee|App.15 A and App. forr us to smoothly move from Nyström-Sinkhom Sinkhorn to sparse Sinkhorn by varying the number of neighbors and landmarks. We can thus freely choose the optimal "operating point" based on the underlying problem and regularization parameter. We discuss the limitations of LCN-Sinkhorn in App. B
for details). It allows
4 THEORETICAL ANALYSIS
Approximation error. The main question we aim to answer in our theoretical analysis is what improvements to expect from adding sparse corrections to Nyström Sinkhorn. To do so, we first analyse approximations of K in a uniform and a clustered data model. In these we use Nyström and LSH schemes that largely resemble k-means, as used in most of our experiments. Relevant proofs and notes for this section can be found in App. Q
Theorem 1. Let Xp and Xq have n samples that are uniformly distributed in d-dimensional closed, locally Euclidean manifold with unit volume. Let furthermore C.p = =laapp Fujlle and Cgg//, Let the landmarks be arranged optimally and regularly so that the expected L2 e = distance to the closest landmark is minimized. Denote R = 1 L.E, || yll2- Assume that the sparse correction K.). = Kij if and only if aaj is one of the nearest neighbors of Epi and that the distance to Tpi k-nearest neighbor dy E R. Then the expected maximum error in row of the LCN approximation KLCN is
BiKics with denoting the index of T 's k-nearest neighbor. Using the upper incomplete Gamma function T(,.-) we can furthermore bound the second term by 2d(T(d) - T(d,2R/M) + +(e (V--3//A/))
=lee
(6)
e -vRRAA S B K S (2R/A) e 222//A)
The error in Eq. (6) is dominated by the first term since 6y R. Note that R only decreases slowly with the number of landmarks since R Z ((dr2)))1/4 2V² Cohn 2017). Moving from pure Nyström to LCN by correcting the nearest neighbors' entries thus provides significant benefits, even for uniform data. For example, by just correcting the first neighbor obtain 68 % improvement in the first term (d (2=32, A= -.005, n= 1000). This is even more pronounced in clustered data.
Theorem 2. Let Xp. X2 S Rdd be distributed inside the same clusters with cluster centers Te. Let be the maximum L2 distance of point to its cluster center and D the minimum distance between two points from different clusters, with 4, D. Let each LSH bucket used for the sparse approximation KPP cover at least one cluster. Let KNys use 1 < 1 S d and KLCN use = 1 optimally distributed landmarks per cluster. Then the maximum error is
le 1e-2/r24 /A
max K- -KNyslloc =
max AE(U.)]
1 1+ (1 -1-1)e--A/ -0c--Di))
(8)
max |K- - K -D/A
(9)
max ||
K- KLCNlla 20 c D/A(1 ( 2r/2 2T/A + D(-)))))
(10)
Since we can lower bound Eq. (8) by 1 - -Dii) we can conclude that the error in K is close to 1 for any reasonably large &. (which is the maximum error possible). The errors in KPP and KLCN on the other hand are vanishingly small, since r << D.
Moreover, these maximum approximation error improvements directly translate improvements in the Sinkhorn approximation. We can show this by slightly adapting the error bounds for an approximate Sinkhorn transport plan and distance due to Altschuler al. (2019).
Theorem 3 (Altcchuler et al. (2019). Let Xp S Idd have n samples. Denote p as the maximum distance between two samples. Let K be an approximation of the similarity matrix K with Kij
[|=====allla?> and K | S -'e-p/ where e' min(1, 506,p log ) When performing
the Sinkhorn algorithm until || pll + 11p'1N all1 S E'2,, the resulting approximate transport plan P and distance d². are bounded by
1d2 d3l S
(KI(PP)P) & 8/A.
(11)
4

---

Under review as a conference paper at ICLR 2021
Convergence rate. We next show that approximate Sinkhorn converges fast regular Sinkhorn by slightly adapting the convergence bound by DDurechensky et al. (2018) to account for sparsity.
Theorem 4 (Dvurechensky et al. (018)). Given the matrix K E 1R and p. the Sinkhorn algorithm gives a transport plan satisfying PPNy - -pl1 + | IN v qll1 S e in iterations
k S 2+ -4 In K Ki > 0f mins, (pin 011
(12)
Backpropagation. Efficient gradient computation almost as important for modern deep learning models as the algorithm itself. These models usually aim at learning the embeddings in Xp and X and therefore need gradients w.r.t. the cost matrix C. We can estimate these either via automatic differentiation of the unrolled Sinkhorn iterations via the analytic solution that assumes exact convergence. Depending on the problem at hand, either the automatic or the analytic estimator will lead to faster overall convergence (Ablin et al.) 2020). LCN-Sinkhorn works flawlessly with automatic backpropagation since it only relies on basic linear algebra (except for choosing Nyström landmarks and L.SH neighbors, for which 20.3)). To enable fast analytic backpropagation we provide analytic gradients in [Proposition 1 Note that both backpropagation methods have runtime linear in the number of points n and m.
simple straight-through estimator (Bengio et al.,
The derivatives of the distances d' and dicx, (Eqs. (1) and (5) and the optimal transport plan P E I w.r.t. the (decomposed) cost matrix C E Ruxm in entropy-regularized OT and LCN-SSikhorn are
Proposition 1.
Od P,
JP.
1 Pipdaeosi X
(13)
OC
0CkI
adieN. aU
dd aw
adhen. θ bog Ksp
Od
= -88
~1(3TU)
= -AP
i log K
-APP
Nys
(14)
PPuM 610d18i1
0 PWwa
a PPus
Pl.ye Pw
Pu. t, P Way
aUk
aUn
OWa
(15)
O W..5. OW
aP:
P Nys.j
PaP
pap
Nys. Su
i log K 39 kl 1j
log K NysKt
s
with dig denoting the Kronecker delta and 1 the Moore-Penrose pseudoinverse. Using these decompo- sitions can backpropagate through LCN-Sinkhorn in time O( m))² + 13).
5 GRAPH TRANSPORT NETWORK
Graph distance learning. The ability to predict similarities or distances between graph- sraph-structured objects useful across a wide range of applications. It can e.g. be used to predict the reaction rate between molecules (Houston et al., 2019), search for similar images (Johnson et al., (015), similar molecules for drug discovery (Birchall et al. 2006), or similar code for vulnerability detection (Li et al. 2019). We propose the graph transport network (GTN) to evaluate approximate Sinkhorn and advance the state the art on this task.
Graph transport network. GTN first uses Siamese graph neural network (GNN) to embed two graphs independently as sets of node embeddings. These embeddings are then matched using enhanced entropy- regropp-regulrrized optimal transport. Given undirected graph G (V,E), with node set V and edge set E, node attributes x; E R H,, and (optional) edge attributes ei, E " with j V, we update the node embeddings in each GNN layer via
h (1) self, a (W(() (1-1) +b()),
node h2
(16)
h
(1) h (4)
sell.i
+ Σ
The selt.
(1
W Wedge
(17)
JEN:
with Nr denoting the neighborhood of node i, h x; h E R Hx for I 2 1, the bilinear layer
W edgee E HyxHH, and the degree normalization mi
(1) and 1i.jj 1/ 1/ /dee, deg, for > 1.
(I)
5

---

Under review as a conference paper at ICLR 2021
This choice of Ta. allows model to handle highly skewed degree distributions while still being able to represent node degrees. We found the choice of non- con-lincarity a not to be critical and chose a LeakyReLU. We do not use the bilinear layer W Ei if there are edge attributes. We aggregate each layer's ' node embeddings to obtain the final embedding of node i
h final
(1)
1 h (2) 11 hn!")
(1)
(18)
Having obtained the embedding sets Hfiimal and Hrfmaall of both graphs we use the L2 distance as a cost function and then calculate the Sinkhorn distance, which is symmetric and permutation invariant w.r.t. the sets Hffinl and Hgmall We obtain the embeddings for matching via h obtain the final prediction via d = di Wout bout, with learnable wou and bow- All weights in GTN are trained and-to-end via backpropagation. For small graphs we use the full Sinkhorn distance and scale to large graphs by leveraging I. LCN--Sinhhorn GTN is more expressive than models that aggegrate node embeddings to a single fixed-size embedding for the entire graph but still scales log-linearly in the number of nodes, as opposed to previous approaches that scale quadratically, Note that GTN inherently performs graph matching and can therefore also be applied to this task.
(h and
Learnable unbalanced OT. Since GTN regularly encounters graphs with disagreeing numbers of nodes it needs to be able to handle cases where ||p|| * ||ql]1 or where not all nodes in one graph have a corresponding node in the other and thus P1m m < p or p"1, < q. Unbalanced OT allows us to handle both of these cases (Peyré & Cuturi 2019). Previous methods did so by swapping these requirements with a uniform divergence loss term on p and q (Frogner al.] 2015 Chizat et al.) 2018). However, these approaches uniformly penalize deviations from balanced OT and therefore cannot adapt to only ignore parts of the distribution. We propose to alleviate this limitation by swapping the cost matrix C' with the bipartite matching (BP) matrix (Riesen & Bunke 2009)
CBp
C
((r.e)
[Cts
i=
cle
cle 0, 000 i itj"
[Ce
i=
(te.g)
C(a,e)
000 i i
(19)
and adaptively computing the costs Ci.e, Cejj and Ca based the input sets Xp and Xq. Using the BP matrix adds minor computational overhead since we only need to save the diagonals Gp.e and Cr of Cp.e and Ce.q. We can then include the additional parts of Cup in the Sinkhorn algorithm [Eq. via
K8pt =
Ki + Cr.e + i Ceaq Ott T
KIPS
KT 8 Cr.q Cr
C 8 m
(20)
where denotes the upper and i the lower part of the vector t. To calculate dž we decompose the transport plan Pse in the same way as CBP, with a single scalar for Pe.s. For GTN we obtain the deletion cost via Ci. = ||o @ æpilla, with a learnable vector & E Rdd
Multi-head OT. Inspired by attention models (Vaswani al. 2017) we further improve GTN by using multiple OT heads. Using K heads means that we calculate late OT in parallel for K separate sets of embeddings representing the same pair of objects and obtain set of distances d R We can then transform these distances to a final distance prediction using a set linear layers h (k) filinal for head A and obtain the final prediction via d = ML.((d). Note that both learnable unbalanced OT and head OT might be of independent interest.
6 RELATED WORK
Log-linear optimal transport. For an overview of optimal transport and its foundations see Peyré Cuturi (2019). On low-dimensional grids and surfaces OT can be solved using dynamical OT (Pa- et [olooooo 2014), convolutions (Solomon et al. (013).. or embedding/hashing
padakis et al. [0044] Solomon
schemes (Indyk& & Thaper 2003] Andoni et al. 2008). In higher dimensions we can use tree- based algorithms (Backurs et al. 2020) hashing schemes (Charikar) 2002). which are however limited to a previously fixed set of points Xq. on which only the distributions and change. For sets that change dynamically (e.g. during training) common method of achieving fogllinear runtime is multiscale approximation of entropy hatropy-regularizer OT (Schmitzer, 2019, Gerber & Maggioni| (017)) Tenetov et al. (2018) recently proposed using a low-rank approximation of the Sinkhorn similarity

---

Under review as a conference paper at ICLR 2021
Table Mcan and standard deviation (w.r.t. last digits, in parentheses) of relative Sinkhorn distance error, IoU of top 0.1 % and correlation coefficient (PCC) of OT plan entries across 5 runs. Sparse Sinkhorn and LCN-Sinkhonn Sinkhorn consistently achieve the best approximation in all 3 measures.
FN-DE
EN-FS
EN-FR
EN-RU
Rel. en. o2
PCC
IoU
Rol. ET. d2
PCC
16U
Rel. erz. <2
PCC
100
Rel. cT.
RCC
160
Factored OT
0.518(1)
0.041:11
0.019/2)
0,552(1)
0
0.02615)
0.336:2)
6.088(1)
0.034:::
0.281(2)
0,055(1) 0.025(2)
Multivcale (163((11)
0.308(14)
6.123(5)
0.645(14)
0.125(12)
0.600(17)
0.3309)
0.121/7)
0.607(16)
0,281(10) 0.125(0)
Nystrmm Sih..
1.18315)
0
(.04(55)
1.176(18)
0.008(1)
0.04866
1.179(18)
0.02244 (.07033)
1.228(18)
0.81(2) 0.47(6)
Sparse Skh. 0.253(()
0.652/4)
0.302())
4.62914) 0.217(1)
0.102(1)
0.220(1)
0.104/2) 0.608(5)
0.272(2) 0.44618) 0.090(1)
LCN-Sinkhoo 0.408(15)
0.673(12)
07197(7) 0.368(12)
0.736(3)
(.201
0.342:9)
0.725(4) 0.209(3)
0.485(10) 0.623(6) 0.210(4)
matrix obtained via a semidiscrete approximation of the Euclidean distance. Altschuler et al. (2019)
improved upon this approach by using the Nyström method for the approximation. These approaches still struggle with high-dimensional real-world problems, as we will show in See. 7
Sliced Wasserstein distance. Another approach to reduce the computational complexity of optimal transport (without entropy regularization) are sliced Wasserstein distances (Rabin et al. 2011)) However, they require the Le distance as cost function and are either unstable in convergence or prohibitively expensive for high- -igh-dimensooal problems (O(nd") (Meng et al. 2019).
Fast Sinkhorn. Another line of work has been pursuing accelerating antropy-regularty OT without changing its computational complexity w.r.t. the number of points. Original Sinkhorn requires O(L/E²) iterations (Dvurechensky et al., 2018) and Jambulapati et al. (2019) recently proposed an algorithm that reduces them to O(1/E). Alaya et al. 2019) proposed to reduce the size of the Sinkhorn problem by screening out neglectable components, which allows for approximation guarantees.
Genevay et al. (2016) proposed using a stochastic optimization scheme instead of Sinkhorn iterations.
Essid & Solomon (2018) and Blondel et al. (2018) proposed alternative regularizations to obtain OT problems with similar runtimes as the Sinkhorn algorithm. This work is largely orthogonal to ours.
Embedding alignment. For overview of cross- cross-lingus lingual word embedding models see Ruder et al. (019)). Unsupervised word embedding aligmment was proposed by Conneau et al. (2018), with subsequent advances by Alvarez-Melis & Jaakkola (018): Grave et al. (2019) Joulin et al. (2018).
Graph matching and distance learning. Most recent approaches for graph matching and graph distance learning either rely on single fixed-dimensional graph embedding (Bai et al.) 2019] Li et. al.] (019). or only use attention or some other strongly simplified variant of optimal transport (Bai et al. 2019 Riba et al., 2018 Li al.] 2019). Others break permutation invariance and are thus ill-suited for this task (Ktena et al.) 20117 Bai et al., 2018). So far only approaches using a single graph embedding allow faster than quadratic scaling in the number of nodes. Compared the Sinkhorn- -baed image model concurrently proposed by Wang et al. (2019) GTN uses no CNN or cross-grapp graph attention, but an enhanced GNN and embedding aggregation scheme. OT has recently been proposed for graph kernels (Maretic et al. [019 Vayer et al. 2019), which can (to limited extent) be used for graph matching, but not for distance learning.
7 EXPERIMENTS
Approximating Sinkhorn. We start by directly investigating different Sinkhorn approximations. To do so we compute entropy-regularized OT on pairs of 000 word embeddings from Conneau al. (2018), which we preprocess with Wasserstein Procrustes alignment in order obtain both close and distant neighbors. We let every method use the same total number of 40 neighbors and landmarks (LCN uses 20 each) and set X == 0.05 (as in Grave et al. (2019)). ). We measure transport plan approximation quality by (a) calculating the Pearson correlation coefficient (PCC) between all entries in the approximated plan and the true P and (b) comparing the sets of 0.1 % largest entries in the approximated and true P using the Jaccard similarity (intersection over union, IoU). In all figures the error bars denote standard deviation across 5 runs. which is often too small to be visible.
Table 1 shows that both sparse Sinkhorn, LCN-Sinkhom and factored OT (Forrow et al.) 2019) obtain distances that are significantly closer to the true d² than Multiscale OT and Nyström-Sinkhan Furthermore, the transport plan computed by sparse Sinkhorn and LCN-Sinkhone show both PCC and IoU that are around twice as high as Multiscale OT, while Nyström-Sinie Sinkhorn and factored OT exhibit almost correlation. CN-Sinkhorn performs especially well in this regard. This is also
7

---

Under review as a conference paper at ICLR 2021
1.0
1.0
1.0
Fact. OT Multse. OT
Sp. Skh. LCN-Skh.
0.8
0.8
0.8
Nys. Skh.
PCC
0.6
PCC
PCC
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0.0
100
200
0.0
100
200
0.0
10
10
10-
10
Runtime (ms)
Neighbors + landmarks
Figure 2: Tradeoff between OT Figure 3: Tradeoff between OT Figure 4: OT plan approximation plan approximation (via PCC) plan approximation and number quality for varying entropy reg- and runtime. Sparse Sinkhorn of- of neighbors/landmarks. LCN-
ularization X. Sparse Sinkhorn performs best for low and LCN- Sinkhorn for moderate and fac- tored OT for very high X.
fers the best tradeoff, with LCN- Sinkhorn achieves the best ap- Sinkhorn trailing closely behind. proximation for low and sparse The arrow indicates factored OT Sinkhorn for high budgets. results far outside the range.
Table 2: Accuracy and standard deviation (w.rt. last digits, in parentheses) across 5 runs for unsuper- vised word embedding alignment with Wasserstein Procrustes. I.CN-Sinkhorn improves upon the original by 3.1 pp. before and 2.0 pp. after iterative CSLS refinement. *Migrated re-run on GPU via PyTorh
Time (s)
EN-ES
ES-EN
EN-FR
FR-EN
EN-DE
DE-EN
EN-RU
RU-EN
Avg.
Original
268
79.2(2)
78.8(2.8)
81.(3)
79.4(9)
7L.72)
66.73311
36.((11)
51.1(111
67.9
Full Sinkhorn
402
81.1()
82.0()
81.20)
81.3()
74.1()
70.7()
37.30)
53.5()
70.1
Multiscale OT
88.2
23.6(31.4)
74.7(3.3)
26.9(31.7) 6.((4.)
35.8(10.4)
47.0(20.5)
0.0()
0.2(1)
26.8
Nyström Skh.
102
64.4(1.0)
59.3(1.2)
64.((1.6) 56.8(4.0)
54.1(6)
47.1(3.5)
14.1(1.2)
22.5(2.4)
47.8
Sparse Skh.
49.2
80.2(2)
81.7(4)
80.9(3)
80.1(2)
72.1(6)
65.1(1.7)
35.5(6)
51.5(4)
68.4
Lin-Sinnh
86.8
81.8(2)
81.3(1.8)
82.0(4)
82.1
73.6(2)
71.3(9)
41.0(8)
55.1(1.4)
71.0
Original ref.
208
81
83.0(3)
82.0(2.5)
83.8(1)
33.014)
77.3(3)
69.7(4.3)
46.00.))
64.0(1.1
72.4
LCN-Skh. + ref.
86.8
81
83.5(2)
83.1(1.3)
83.8(2)
83.6(1)
77.2(3)
72.8(7)
51.((2.)
59.2(1.9)
74.4
evident in Fig. which shows how the 10a x 10a approximated OT plan entries compared to the true Sinkhorn values.
Fig. 2 shows that sparse Sinkhorn offers the best trade-off between runtime and OT plan qual- ity. Factored OT exhibits a runtime to 10 times longer than the competition due its iterative refinement scheme. LCN-Sinkhorn performs best for use cases with constrained memory (few neigh- bors/landmarks), as shown Fig. 3 The number of neighbors and landmarks directly determines memory usage and is linearly proportional to the runtime (see App. K). Fig. shows that sparse Sinkhorn performs best for low regularizations, where LCN-Sinkhorn fails due the Nyström part going out of bounds. Nyström Sinkhorn performs best at high values and LCN-Sinkhorn always performs better than both (as long as it can calculated). Interestingly, all approximations except factored OT seem to fail at high X. We defer analogously discussing the distance approximation to and dataset
App. L. All approximations scale linearly both in the number of neighbors/landmall size, shown in App. K. Overall, we see that sparse Sinkhorn and LCN- Sinkhorn yield significant improvements over previous approximations. However, do these improvements also translate to better performance on downstream tasks?
Embedding alignment. Embedding alignment is the task of finding the orthogonal matrix R E Rddd that best aligns the vectors from two different embedding spaces, which is e.g. useful for unsupervised word translation. We use the experimental setup established by [onneau et al. (2018) by migrating Grave al. (2019)'s 's implementation to PyTorch. The only change we make is using the full set of 20000 word embeddings and training for 300 steps, while reducing the learning rate by half every 100 steps. We do not change any other hyperparameters and do not use unbalanced OT. After training we match pairs via cross-domain similarity local scaling (CSLS) (Conneau et al. 2018). We use 10 Sinkhorn iterations, 40 neighbors for sparse Sinkhorn, and 20 neighbors and landmarks for LCN-Sinkhorn (for details see/App. H). We allow both multiscale OT and Nyström Sinkhorn to as many landmarks and neighbors as can fit into GPU memory and finetune both methods.
Table shows that using full Sinkhorn yields significant improvement in accuracy on this task
8

---

Under review as a conference paper at ICLR 2021
Table 3: RMSE for GED regression across
Table 4: RMSE for graph distance regression across
3 runs and the targets' standard deviation a.
3 runs. Using LCN-Sinkhorn with GTN increases the error by only 10 % and allows log- log-linear scaling.
GTN outperforms previous models by 48 %.
Linux
AIDS30
Pref. att.
GED
PM [10-2]
0.184
16.2
48.3
AIDS30
Pref. att. Pref. att. 200
SiamMPNN
0.090(7)
13.8(3)
12.1(6)
16.2
48.3
10.2
SimGNN
0.039
4.5(3)
8.3(1.4)
Full Sinkhorn
3.7(1)
4.5(3)
1.27(6)
GMN
0.015()
10.3(6)
7.8(3)
Nyström Skh.
3.6(3)
6.2(6)
2.43(7)
GTN, head 0.022(1)
3.7(1)
4.5(3)
Multiscale OT
11.2(3)
27.4(5.4)
6.71(44)
8 OT heads 0.012(1)
3.2(1)
3.6(2)
Sparse Skh.
44.0(30.4)
40.7(8.1)
7.57(1.09)
Balanced OT 0.034(1)
15.3(1)
27.4(9)
LCN-Skh.
4.0(1)
5.1(4)
1.4((15)
compared to the original approach of performing Sinkhorn on randomly sampled subsets of embed- dings (Grave et al. 2019). LCN-Sinkhorn even outperforms the full version in most cases, which is likely due to regularization effects from the approximation. It also runs 4.6x faster than full Sinkhorn and 3.1x faster than the original scheme. Sparse Sinkhorn runs 1.8x faster than LCN-Sinkhorn but cannot match its accuracy. LCN-Sinkhorn still outcompetes the original method after refining the embeddings with iterative local CSLS (Conneau et al. 2018). Both multiscale OT and Nyström Sinkhorn fail at this task, despite their larger computational budget. This shows that the improvements achieved by sparse Sinkhorn and LCN-Sinkhorn have an even larger impact in practice.
Graph distance regression. The graph edit distance (GED) is useful for various tasks, such as image retrieval (Xiao al. 2008) or fingerprint matching (Neuhaus & Bunke 2004), but its computation is NP-complete (Bunke & Shearer 1998). Therefore, to use it on larger graphs we need to learn an approximation. We use the Linux dataset by Bai et al. (2019) and generate new datasets by computing the exact GED using the method by| Lerouge et al. (2017) on small graphs (s 30 nodes) from the AIDS dataset (Riesen & Bunke 2008) and set of preferential attachment graphs. We compare GTN to 3 state-of-the-art baselines: SiameseMPNN (Riba et al. 2018). SimGNN (Bai et al. 2019), and the Graph Matching Network (GMN) (Li et al. 2019). We tune the hyperparameters of all baselines and GTN the validation set via a grid search. For more details see App. H| We first test both GTN and the proposed OT enhancements, Table 3 shows that GTN improves upon competing models by 20 % with a single head and by 48 % with 8 OT heads. These improvements break down when using regular balanced OT, showing the importance of learnable unbalanced OT. Having established GTN as a state-of-the-art model formance when using approximate OT. To test this we additionally generate a set of larger graphs with around 200 nodes and use the Pyramid matching (PM) kernel (Nikolentzos et al. 2017) as the prediction target, since these graphs are too large to compute the GED. See App. J for hyperparameter details. Table sjshows that both sparse Sinkhorn and the multiscale method using 4 (expected) neigh- bors fail at this task, demonstrating that the low-rank approximation in LCN has a crucial stabilizing effect during training. Nyström Sinkhorn with 4 landmarks performs surprisingly well on the AIDS30 dataset, suggesting
next ask whether we can sustain its per-
overall low-rank structure with Nyström acting as regularization. However, it does not perform as well on the other two datasets. Using LCN- CCN--inkhorn with 2 neighbors and landmarks works well on all three datasets, with RMSE increased by only 10 % compared to full indeed scales linearly in the number
GTN. App. Kffurthemoore shows that GTN with of nodes across multiple orders of magnitude. This model thus allows perform graph matching and distance learning on graphs that are considered large even for simple node-level tasks (20 000 nodes).
8 CONCLUSION
Locality sensitive hashing (LSH) and the novel locally corrected Nyström (LCN) method enable fast and accurate approximations entropy-regullrizzed OT with lon-linear runtime: Sparse Sinkhorn and LCN-Sinkhorn. The graph transport network (GTN) is one example for such model, which can be substantially improved with learnable unbalanced OT and multi-head OT. It sets the new state of the art for graph distance learning while still scaling log- linearlyll with graph size. These contributions enable new applications and models that are both faster and more accurate, since they can sidestep workarounds such as pooling.

---

Under review as a conference paper at ICLR 2021
REFERENCES
Pierre Ablin, Gabriel Peyré, and Thomas Moreau, Super-efficieney of automatic differentiation for functions defined as a minimum. In ICML, 2020.
Mokhtar Z. Alaya, Maxime Berar, Gilles Gasso, and Alain Rakotomamonjy. Screening Sinkhorn Algorithm for Regularized Optimal Transport. In NeurIPS, 2019.
Jason Altschuler, Francis Bach, Alessandro Rudi, and Jonathan Niles-Weeks Massively scalable Sinkhorn distances via the Nyström method. In NeurIPS, 2019.
David Alvarez-Melis and Tommi S. Jaakkola. Gromov-Wasserstein Alignment of Word Embedding Spaces. In EMNLP, 2018.
Alexandr Andoni, Piotr Indyk, and Robert Krauthgamer. Earth mover distance over high-dimensional spaces. In ACM-SIAM symposium on Discrete algorithms (SODA), 2008.
Alexandr Andoni, Piotr Indyk, Thijs Laarhoven, Ilya P. Razenshteyn, and Ludwig Schmidt. Practical and Optimal LSH for Angular Distance. In NeurIPS, 2015.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein Generative Adversarial Networks. In ICML, 2017.
Arturs Backurs, Yihe Dong, Piotr Indyk, Ilya Razenshteyn, and Tal Wagner. Scalable Nearest Neighbor Search for Optimal Transport. In ICML,, 2020.
Yunsheng Bai, Hao Ding, Yizhou Sun, and Wei Wang. Convolutional Set Matching for Graph Similarity. In NeurIPS-W, 2018.
Yunsheng Bai, Hao Ding, Song Bian, Ting Chen, Yizhou Sun, and Wei Wang. SimGNN: A Neural Network Approach to Fast Graph Similarity Computation. In WSDM, 2019.
Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating or Propagating Gradients Through Stochastic Neurons for Conditional Computation. CoRR, 1308.3432. 2013.
Christian Berg, Jens Peter Reus Christensen, and Paul Ressel. Harmonic Analysis on Semigroups. Number 100 in Graduate Texts in Mathematics. 1984.
Kristian Birchall, Valerie Gillet, Gavin Harper, and Stephen D. Pickett. Training Similarity Measures for Specific Activities: Application to Reduced Graphs. Journal of Chemical Information and Modeling, 46(2):577-586, 2006.
Mathieu Blondel, Vivien Seguy, and Antoine Rolet. Smooth and Sparse Optimal Transport. In AISTATS, 2018.
Olivier Bousquet, Sylvain Gelly, llya Tolstikhin, Carl-Johann Simon-Gabriel, and Bernhard Schoelkopf. From optimal transport to generative modeling: the VEGAN cookbook. arXiv, 1705.07642, 2017.
Horst Bunke and Kim Shearer. A graph distance metric based on the maximal common subgraph. Pattern Recognition Letters, 19(3):55-2255 1998.
Moses Charikar. Similarity estimation techniques from rounding algorithms. In ACM symposium on Theory of computing (STOC), 2002.
Lénaïc Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard. Scaling algorithms for unbalanced optimal transport problems. Mathematics of Computation, 87(114)2263-2609, 2018.
Henry Cohn. A Conceptual Breakthrough in Sphere Packing. Notices of the American Mathematical Society, 64(02):102-115, 2017.
Alexis Conneau, Guillaume Lample, Mare' Aurelio Ranzato, Ludovic Denoyer, and Hervé Jégou. Word translation without parallel data. In ICLR, 2018.
10

---

Under review as a conference paper at ICLR 2021
Nicolas Courty, Rémi Flamary, Amaury Habrard, and Alain Rakotomamonjy. Joint distribution optimal transportation for domain adaptation. In NeurIPS, 2017.
Marco Cuturi. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. In NeurlPS, 2013.
Pavel E. Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational Optimal Trans- port: Complexity by Accelerated Gradient Descent Better Than by Sinkhorn's 's Algorithm. In ICML, 2018.
Montacer Essid and Justin Solomon. Quadratically Regularized Optimal Transport on Graphs. SIAM Journal on Scientific Computing, 40(4):A1961-A195)) 2018.
Matthias Fey and Jan E. Lenssen. Fast Graph Representation Learning with PyTTrch Geometric. In ICLR-W, W, 2019.
Aden Forrow, Jan- Christian Hütter, Mor Nitzan, Philippe Rigollet, Geoffrey Schiebinger, and Jonathan Weed. Statistical Optimal Transport via Factored Couplings. In AISTATS, 2019.
Charlie Frogner, Chiyuan Zhang, Hossein Mobahi, Mauricio Araya-Polo, and Tomaso A. Poggio. Learning with Wasserstein Loss. In NewrIPS, 2015.
Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis R. Bach. Stochastic Optimization for Large-scale Optimal Transport. NeurlPS, 2016.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning Generative Models with Sinkhorn Divergences. In AISTATS, 2018.
Samuel Gerber and Mauro Maggioni. Multiscale Strategies for Computing Optimal Transport Mach. Learn. Res., 18:72:1-72:32, 2017.
Edouard Grave, Armand Joulin, and Quentin Berthet. Unsupervised Alignment of Embeddings with Wasserstein Procrustes. In AISTATS, 2019.
Paul L. Houston, Apurba Nandi, and Joel M. Bowman. A Machine Learning Approach for Prediction of Rate Constants. The Journal of Physical Chemistry Letters, 10(17)52250-528, 2019.
Piotr Indyk and Nitin Thaper. Fast image retrieval via embeddings. In ICCV-W, 2003.
Arun Jambulapati, Aaron Sidford, and Kevin Tian. A Direct tilde (O)((eesilll Iteration Parallel Algorithm for Optimal Transport. In NeurlPS, 2019.
Justin Johnson, Ranjay Krishna, Michael Stark, Li-Jia Li, David A. Shamma, Michael S. Bernstein, and Fei-Fei Li. Image retrieval using scene graphs. In CVPR, 2015.
Armand Joulin, Piotr Bojanowski, Tomas Mikolov, Hervé Jégou, and Edouard Grave. Loss in Translation: Learning Bilingual Word Mapping with Retrieval Criterion. EMNLP, 2018.
Soña Ira Ktena, Sarah Parisot, Enzo Ferrante, Martin Rajchi, Matthew Lee, Ben Glocker, and Daniel Rueckert. Distance Metric Learning Using Graph Convolutional Networks: Application to Functional Brain Networks. In MICCAI, 2017.
Julien Lerouge, Zeina Abu- Absheh, Romain Raveaux, Pierre Héroux, and Sébastien Adam. New binary linear programming formulation to compute the graph edit distance. Pattern Recognit. 72: 254 265, 2017.
Yujia Li, Chenjie Gu, Thomas Dullien, Oriol Vinyals, and Pushmeet Kohli. Graph Matching Networks for Learning the Similarity of Graph Structured Objects. In ICML, 2019.
Hermina Petric Maretic, Mircille El Gheche, Giovanni Chierchia, and Pascal Frossard. GOT: An Optimal Transport framework for Graph comparison. In NewrlPS, 2019.
Cheng Meng. Yuan Ke, Jingyi Zhang, Mengrui Zhang, Wenxuan Zhong, and Ping Ma. Large-scale optimal transport map estimation using projection pursuit. In NewrIPS, 2019.
11

---

Under review as a conference paper at ICLR 2021
Cameron Musco and Christopher Musco. Recursive Sampling for the Nystrom Method. In NewrIPS, 2017.
Michel Neuhaus and Horst Bunke. An Enror-Tolerant Approximate Matching Algorithm for Attributed Planar Graphs and Its Application to Fingerprint Classification. In Structural, Syntactic, and Statistical Pattern Recognition, 2004.
Giannis Nikolentzos, Polykarpos Meladianos, and Michalis Vazirgiannis. Matching Node Embed- dings for Graph Similarity. In AAAI, 2017.
David Nistér and Henrik Stewénius. Scalable Recognition with Vocabulary Tree. In CVPR, 2006.
Nicolas Papadakis, Gabriel Peyré, and Édouard Oudet. Optimal Transport with Proximal Splitting. SIAM J. Imaging Sciences, (1):2212-238, 2014.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang. Zachary Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Perffrmance Deep Learning Library. In NeurlPS, 2019.
Loic Paulevé, Hervé Jégou, and Laurent Amsaleg. Locality sensitive hashing: comparison hash function types and querying mechanisms. Pullern Recognit. Lett., 31(1 11:1138-113355 2010.
Allon G Percus and Olivier C Martin. Scaling Universalities ofkth-Nearest Neighbor Distances on Closed Manifolds. Advances in Applied Mathematics, 21(3):424-45) 436, 1998.
Gabriel Peyré and Marco Cuturi.
Computational Optimal Transport. Foundations and Trends in Machine Learning, 1(5-6))355-677, 2019.
Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein Barycenter and Its Application to Texture Mixing. In Scale Space and Variational Methods in Computer Vision (SSVM), 2011.
Pau Riba, Andreas Fischer, Josep Lladós, and Alicia Fornés. Learning Graph Distances with Message Passing Neural Networks. In ICPR, 2018.
Kaspar Riesen and Horst Bunke. IAM Graph Database Repository for Graph Based Pattern Recogni- tion and Machine Learning. In Structural, Syntactic, and Statistical Pattern Recognition, 2008.
Kaspar Riesen and Horst Bunke. Approximate graph edit distance computation by means of bipartite graph matching. Image Vis. Comput., 27(7):950-9)) 959, 2009.
Sebastian Ruder, Ivan Vulic, and Anders Søgaard. A Survey of Cross-lingual Word Embedding Models. J. Artif. Intell. Res., 65:569-631, 2019.
Bernhard Schmitzer. Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems. SIAM Journal Scientific Computing, (143):A1443-A1444 2019.
Anshumali Shrivastava and Ping Li. Asymmetric LSH (ALSH) for Sublinear Time Maximum Inner Product Search (MIPS). In NeurIPS, 2014.
Justin Solomon, Raif M. Rustamov, Leonidas J. Guibas, and Adrian Butscher. Earth mover's distances on discrete surfaces. ACM Trans. Graph., 33(4)667:1-67112, 2014.
Justin Solomon, Fernando de Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas J. Guibas. Convolutional wasserstein distances: efficient optimal transportation on geometric domains. ACM Trans. Graph., 34(4):66: 1-66: 2015.
Robert E. Tarjan. Dynamic trees as search trees via euler tours, applied to the network simplex algorithm. Mathematical Programming, 1892) 1997.
Evgeny Tenetov, Gershon Wolansky, and Ron Kimmel. Fast Entropic Regularized Optimal Transport
Using Semidiscrete Cost Approximation. SIAM J. Sei. Comput,, 4065):A3400-A3422 2018.
12

---

Under review as a conference paper at ICLR 2021
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention All you Need. In NeurIPS, 2017.
Titouan Vayer, Nicolas Courty, Romain Tavenard, Laetitia Chapel, and Rémi Flamary. Optimal Transport for structured data with application graphs. In ICML, 2019.
Jingdong Wang. Heng Tao Shen, Jingkuan Song, and Jianqiu Ji. Hashing for Similarity Search: Survey. CoRR, 1408.2927, 2014.
Runzhong Wang. Junchi Yan, and Xiaokang Yang. Learning Combinatorial Embedding Networks for Deep Graph Matching. In ICCV, 2019.
Christopher K. I. Williams and Matthias Seeger. Using the Nyström Method to Speed Up Kernel Machines. In NeurIPS, 2001.
Bing Xiao, Xinbo Gao, Dacheng Tao, and Xuelong Li. HMM-based graph edit distance for image indexing. Int. J. Imaging Systems and Technology, 18(2-3):209-22 -200 2008.
Kai Zhang, Ivor W. Tsang, and James T. Kwok. Improved Nyström low-rank approximation and error analysis. ICML, 2008.
13

---

Under review as a conference paper at ICLR 2021
A COMPLEXITY ANALYSIS
Sparse Sinkhorn. common way of achieving a high p1 and low P2 in L.H is via the AND-OR construction. In this scheme we calculate B.r hash functions, divided into B sets (hash bands) of hash functions each. A pair of points is considered as neighbors if any hash band matches completely. Calculating the hash buckets for all points with b hash buckets per function scales as O((n+ ((n+m) m Bbr for the hash functions we consider. As expected, for the tasks and hash functions we investigated we obtain approximately m/ m/6" and n/b" neighbors, with b" hash buckets per band. Using this we can fix the number of neighbors to a small, constant 8 in expectation with b" = min(n, m3.3. We thus obtain sparse cost matrix CPP with O (max(n, m)) non-infinite values and can calculate and in linear time (i max(n, m))3,, where Niink S 2 + m (K,yl K220) mini, (pr.g,)) (see Theorem 4) denotes the number of Sinkhorn iterations. Calculating the hash buckets with r log min[n, mlo takes ( + m) m)dBBbllog min(n, m) - log 8)/ log b). Since B, b, and | are small, we obtain roughly log-linear scaling with the number of points overall, i.e. O(n log for n R m..
log
LCN-Sinkhorn. Both choosing landmarks via k-means++ sampling and via k-means with fixed number of iterations have the same runtime complexity of ((( (m)mm))d) m)ld). Precomputing W can be done in time Ond (s). The low-rank part of updating the vectors 8 and t can be computed in Onn + 12 + Im), with chosen constant, i.c. independently of n and m. Since sparse Sinkhorn with LSH has a log-linear runtime we again obtain log-lineat overall runtime for CCN-Sinkhorn.
B LIMITATIONS
Sparse Sinkhorn. Using a sparse approximation for K works well in the common case when the regularization parameter À is low and the cost function varies enough between data pairs, such that the transport plan P resembles sparse matrix. However, it can fail if the cost between pairs is very similar or the regularization very high, if the dataset contains many hubs, i.c. points with large number of neighbors, or if the distributions p q are spread very unevenly. Furthermore, sparse Sinkhorn can be too unstable to train a model from scratch, since randomly initialized embeddings largely alleviates these limitations.
often have no close neighbors (see Sec. 7)
LCN-Sinkhorn. Since we cannot calculate the full cost matrix, LCN-Sinkhorn cannot provide accuracy guarantees in general. Highly concentrated distributions p and q might have adverse effects on LCN-Sinkhorn. However, we can compensate for these by sampling landmarks or neighbors proportional to each point's probability mass. We therefore do not expect down in this scenario. If the regularization parameter is low or the cost function varies greatly, we sometimes observed stability issues (over- and underflows) with the Nyström approximation because of the inverse A 11 which cannot be calculated in log-space. Due to its linearity the Nyström method furthermore sometimes approximates similarities as negative values, which leads to a failure if the result of the matrix product in Eq. (2) becomes negative, In these extreme cases we also observed catastrophic elimination caused by the correction KT- Since this essentially means that optimal transport will be very local, we recommend using sparse Sinkhorn in these scenarios. This again demonstrates the complementarity of the sparse approximation and Nyström: In cases where one fails we can often resort to the other.
to break
C PROOF OF [HEOREM
We first prove a lemma that will be useful later on.
Lemma A. Let K be the Nyström approximation of the similarity matrix Kij =ll Let I; and xy be data points with equal L2 distance F; and ry to all landmarks, which have the same distance A to each other. Then
le (retrs///
KG 1+ (1 10/
(21)
Proof. The inter-landmark distance matrix
A 1/1 + (1 )/
(22)
14

---

Under review as a conference paper at ICLR 2021
where 1xx denotes the constant 1 matrix. Using the identity
bb
1
(bluxn + (a - b) (a - (a + (n 1)6) 1axn +
a- In
(23)
we compute
K2 - U,AAAVV V.s
/-/s/>>
= (c E ...) (1- + (1 1)e- lixt + A-e-A/x A/X 11
--A
e
ennttry)/
1+ (1 - 1)e A/A +1) = 1 1-- A/A 1 (1 - 1)e 2/7
-1/
((rttr)/ 1 1-
1-e-A/A A/A
le (1---/1
1+ (1-1)e-22 1)e-A/A
(24)
Now consider the error K,, KLCN,i,lla The k - nearest neighbors are covered by the sparse correction and therefore the next nearest neighbor has distance óg. The expected distance from the closest landmark is greater than the expected distance inside the surrounding d-ball of radius R, i.e. E[r] 2 Ev(n) ir] = d+1 R. Because furthermore 6µ < R. the error is dominated by the first term and the maximum error in row given by the k-nearest neighbor i, denoted by j. Thus
[|[| K;,- :I..,i.llss] R Kii KI R k.g] E, i. R(e 6e/A] IE| LLN
(25)
Note that we can lower bound the first term using Jensen's inequality. However, we were unable to find a reasonably tight upper bound and the resulting integral (ignoring exponentially small boundary effects, see Pereus & Martin (1988)
Gd/23)1/4
n!
dV dr
Ele 66-
(1-
dr,
(26)
(n k)( - 1)!
with the volume of the d-ball
4/2
V(r) (d/2)!
(27)
does not have an analytical solution. We thus have to resort to calculating this expectation numerically.
We lower bound the second term by (1) ignoring every landmark except the closest one, since additional landmarks can only increase the estimate KLCN,i.j. We then (2) upper bound the L2 distance to the closest landmark by vdR/2, since this would be the furthest distance to the closest point in a d-dimensional grid. Any optimal arrangement minimizing Efmingez yez || ylla x E X,] would be at least as good as a grid and thus have furthest distances as small or smaller than those in grid. Thus,
(1) E E 2 2r/2 ( Van/A
(28)
We upper bound this expectation by considering that any point outside the inscribed sphere of the space closest to a landmark (which has radius Rt would be further away from the landmarks and thus have a lower value e d// We can therefore reduce the space over which the expectation is taken to the ball with radius R,, i.c.
E KiK.N, S By(R) (i.N.i.)
(29)
Next we (1) ignore the contributions of all landmarks except for the closest 2, since a third landmark must be further away from the data point than v3R, adding error of (c We then (2)
15

---

Under review as a conference paper at ICLR 2021
lower bound the distances of both points to both landmarks by the closest distance landmark
r = min [ - 21, 21 ||
212ll2.
||,
21, 12, ||, T12 2 and use Lemma to obtain
Ev(R) (L
(V(R) [KLNN landmarks, i, + O 2V/3/A)
(E--2//RAA)
(2)
2e 2e/2r/A
(v(
1 1+e 1+e-2R1 2R/A
+ + O -2/3//A/
(30)
2Ev(m)
22/211
O (V/3RAA)) 1 2R/A
Assuming Euclideanness (I)) we obtain
Evenyle le-0l2-2121- = V(R) L. e 2r/A dv(r) dr dr = d R
R / e e-2r/2,d-100 dr
{"
(31)
(2R/Ajd (L(d) r(d, 2R/
D PROOF OF THEOREM 2
Note that this theorem does not use probabilistic arguments but rather geometrically analyzes the maximum possible error. KTP is correct for all pairs inside cluster and otherwise. We therefore obtain the maximum error by considering the closest possible pair between clusters. By definition, this pair has distance D and thus
max K
K e-DAA
(32)
LCN is also correct for all pairs inside cluster, again consider the closest possible pair æ; æj between clusters. We furthermore only consider the landmarks of the two concerned clusters, adding an error of Oee-D/A) Hence,
KLCN, if = (e- (U(DD// (: 1 e-(artDi/A)
(2-11D))/
+/A
- U1+22)/A (e-r/A -(++D)/A)
1
(2r1
D1/2
1/a
(
1 1- (4r+22D (C-A/A c-N++D/A)
r/A _ (3r 12D1/2
1 1 (1r+222/A (3- - (4+1))/A + (2r+1)//> - (144331))
(2+D))/A
1 (1r+22// (2 - 2r/2 _ (22+22)/22
DItE (2/2 - 2//2- -O(e 2D/a)
(33)
and thus
max K - - -2- D(-))))).
(34)
For pure Nyström we need to consider the distances inside a cluster. In the worst case two points overlap, i.e. K,j 1, and lie the boundary of the cluster. Since < D we again only consider the landmarks the concerned cluster, adding an error of O(e O(e-D/A) Because of symmetry we can optimize the worst-case distance from all landmarks by putting them on a (1 ))-smpleex centered on the cluster center. Since there are most landmarks in each cluster there is always one direction in which the worst-case points are away from all landmarks. The circumradius of an (1 1)-simplex
with side length A is / 1 A. Thus, the maximum distance to all landmarks is //2 2 2
Using Lemma A we therefore obtain the Nyström approximation
KNys,8 = le 2 2 /
1+ (1-1)) A/A + + D/A)
(35)
n
16

---

Under review as a conference paper at ICLR 2021
E NOTE ON THEOREM 3
Lemmas C-F and and thus Theorem by Altschuler al. (2019) also valid for Q outside the simplex so long ||Q]| and only has non-negative entries. Any P returned by Sinkhorn fulfills these conditions. Therefore the rounding procedure given by their Algorithm 4 is not necessary for this result.
Furthermore, to be more consistent with Theorems 1 and 2 we use the L2 distance instead of L3 in this theorem, which only changes the dependence p.
F NOTES ON THEOREM
To adapt Theorem by Dvurechensky et al. (2018) to sparse matrices (i.e. matrices with some Kij = 0) need to redefine
V min
K Kj > 0),
(36)
1.
i.e. take the minimum only w.r.t. non-zero elements in their Lemma
G PROOF OF PROPOSITION
Theorem A (Danskin's theorem). Consider a continuous function R x 2 - R. with the
compact set C R. If ff z) is convex x for every Z and Φ(x, 2) has a unique maximizer 2, the derivative of
f(z) max (xx 2)
(37)
is given by the derivative at the maximizer, i.e.
0f dō(x,
(38)
Ōx
We start by deriving the derivatives of the distances. To show that the Sinkhorn distance fulfills the conditions for Danskin's theorem first identify (P))) We next observe that the restrictions P1,, p and PT 1, q define a compact, convex set for Furthermore, & continuous function and linear C, i.e. both convex and concave for any finite P. Finally, dCC,, P) is concave in P since (P. C) is linear and AH(P) is concave. Therefore the maximizer P unique and Danskin's theorem applies the Sinkhorn distance. Using
C, z P, and @CC, P) -
@CNy..5 θ OUN
Wi
Wij
-Alog( Σ Uia Waj)
= -15 Σ. U
(39)
OUNl
Woj
K
0CNys,5
Uik
U
Aloog( Σ Uha Way)
-681
- -68s K
(40)
OWHl
9Wm
Σ. U W...O
E, Pu. Pw..j
Bel, Σ, UN Woj
U.U WoJ
3,t, Σn Uno Wai
3,t,
(41)
E. Uia Waj
En Uio Waj
17

---

Under review as a conference paper at ICLR 2021
and the chain rule we can calculate the derivative w.r.t. the cost matrix
Od 0 ( (- (-(, ((())) =,,
(42)
ac
ac
ad[eN.
Σ dCNys.ij ddicN,e =
8 Wij Piys,jj
Σ
OUkI
OUkI @CNys,t
KNys.is
(43)
=
Σ
8 Wij3,t, = -33k 2 Wists = (->A(WWB))
1
aCNys.i ōdyoN. = - Σ
Pays,d
Uil
OW
O @CNys,ij
KNys.ij
1.
(.
(44)
Σ UnUnd - Σ
8;Uik di 1->~(3-UU)T"))
adfen. and
Bd Flng K Nys
0d.. 8C
and
directly follow from
Flog K
Theorem (Implicit function theorem). Let f Th' Rm R be a continuously differentiable function with f(a, b) O. If its Jacobian matrix J f((ma))) as. (a. b) is invertible, then there ay: exists an open set a E U TR' on which there exists a unique continuously differentiable function g: U + m with g(a) and Vo U ((xxg(()) = Moreover
0g. dx,
( (( (
Ix ( f(x,, ( ((x)((())))))))))))))
(45)
:
Next we derive the transport plan derivatives. To apply the implicit function theorem we identify = C and y = as flattened matrices. We will index these flat matrices via index pairs to simplify interpretation. We furthermore identify
S(C, = Op (P.CCH- = C + X(log P ++1)
(46)
The minimizer P cannot lie on the boundary of P's valid region since limp- ap logp
limp-+0 log, + 1 -co and therefore Pi > 0. Hence, f(C, P) 0 with P(C)
arg mimppppPe (P. C) AH(P) and we find that g(x) = P(C). We furthermore obtain
(C, J.,(C,)),P. (C, (CC,((())))))))))) OP.b (Ciz + A(log P.j. + 1)) = SeadsoA/P.j.
(47)
0CkI (c, aCH
(Gat + X Pat 1)) Sak-bss
(48)
J, P., P(C) is hence diagonal matrix and invertible since A/P.j > 0. We can thus use the implicit function theorem and obtain
aP, (UI..C.P),P aC
(as(C, P(C))
o.b
(49)
Eonadee - f - P,yOs1Op1
-
b.b
To extend this result LCN-OT
aPuM
Σ Pays, Pt. W.kj = 6ia Pt Wi
(50)
0 Piys.
ONNys,ab
i OPwss
Σ
Pt Pryss 65 Ui.in
(51)
OP Nys.ab
PPNys,ab
18

---

Under review as a conference paper at ICLR 2021
and the chain rule obtain
PPuM = Σ
&Cyssedd 0 @ ayy PP
0UkI
aUK
@CNys.od P Nyys,ab
a,b,o,d
-xdek KNys,ca/ )
Wid
(-1 Pryrasdocbbs)) (diePirs)) Pircss)
Σ abbind
(52)
PNys.ib Pi
E
Wis
P WibjOsk 8sk Σ 8; Wist du Σ Pw..oPi..s
Kyysab
8stdg184,
aPw.j
@C a ab aPw.ys 0Un aCNys.od PNys,ab
Σ a.bad,
aU
Wid
)
(-) Prysaasdaddad) (ssoPl.ce)
Σ abbadd
Adck KNys,od/
(53)
Pays.ki pt Pliuk Wist, Wisty Pl Pt Pi Wij
KNyskj
P.u.j
@Cyys,od ONyys,ab OPv.j
Σ
OW Rl
OW aCNys.ed PNyss,bb
a.b.c,d
Uok
) (-) Prsssssbucru)
Σ
-ddai KNys,od/
(54)
abbind
Uik Priy,il Piays s: u Pw.1 Pu,at P Piras Kyys,st
iPw.ss
OCNys.od PNyys,ab OPw.is
Σ
OW 3CNys,ed ò PNys,,ab
abboad
Uck
(-2 PNyaoddacol (oxpPttin)
Σ abb.cdd
-Xddt KNys,od
(55)
Physajj p 8, Σ SaUntt, 6git, Pl,aa Pu Pu.a
Eaa KNys,as
BPP Flog K and
aP Foog K directly follow from ap SC.
We can calculate the pseudoinverses Pjj (PT and (P in time O(n + m) + (3) since Py E Æaxt and Pw Rixm We do not fully instantiate the matrices required for backpropagation but instead save their decompositions, similar to the transport plan PNys P& Pw. We can then compute backpropagation in time O(nn m))²) by applying the sums over and in the right order. We thus obtain O(nn + m))² 13) overall runtime for backpropagation.
H CHOOSING LSH NEIGHBORS AND NYSTRÖM LANDMARKS
We focus on two L.S methods for obtaining near neighbors. Cross-polytope LSH (Andoni al. 2015) uses random projection matrix R dxb/22 with the number of hash buckets b, and then decides on the hash bucket via h(x) = arg max((a* R x R)), where denotes concatenation. K-means LSH computes k- k-means and uses the clusters as hash buckets.
We further improve the sampling probabilities of cross-poly pross-polytope LSH via the AND-OR construction. In this scheme calculate B hash functions, divided into B sets (hash bands) of hash functions
19

---

Under review as a conference paper at ICLR 2021
Table Graph dataset statisties.
Distance (lest set)
Graphs train/val/test
Avg. nodes
Avg. edges per graph
Node types
Edge
Gruph type Molecules
Distance
Meun
Std. dev.
per graph
types
AIDS30
GED
50.5
16.2
144/48/48
21.6
44.6
53
Linux
Program dependence
GED
0.567
0.181
600/200/200
7.6
6.9
Pref, att.
Initial attrattiveness
GED
106.7
48.3
144/48/48
21.6
75.4
Pref, att, 200
Initial uttractiveness
PM
0.400
0.102
144/48/48
199.3
938.8
Pref.
Initial attractiveness
PM
0.359
0.163
144/48/48
2045.6
11330
Pref. att. 20%
Initial attractiveness
PM
0.363
0.151
144/4848
20441
90412
each. A pair of points is considered as neighbors if any hash band matches completely. K-means LSH does not work well with the AND-OR construction since its samples are highly correlated, For large datasets we use hierarchical k-means instead (Paulevé et al., 2010 Nistér & Stewénius (0066.
Since the graph transport network (GTN) uses the L2 distance between embeddings as a cost function we use (hierarchical) k-means LSH and k-means Nyström in both sparse OT and LCN- OT. For embedding alignment we use cross-polytope I.SH for sparse OT since similarities are measured via the dot product. For LCN-OT we found that using k-means LSH works better with Nyström using k-means++ sampling than cross-polytope LSH. This is most likely due to better alignment between LSH samples and Nyström. We convert the cosine similarity to a distance via
doos
(Berg et al. 1984) to use k-means with dot product similarity. Note that 22
this is actually based on cosine similarity, not the dot product. Due to the balanced nature of OT we found this more sensible than maximum inner product search (MIPS). For both experiments we also experimented with uniform and recursive RLS sampling but found that the above mentioned methods work better.
I IMPLEMENTATIONAL DETAILS
Our implementation runs in batches on GPU via PyTorch (Paszke et al. 2019) and PyTorch Scatter (Fey & Lenssen 2019). To avoid over- and underftows use log-stabilizztion throughout, i.e. save all values in log-space and compute all matrix- vector products and additions via the log-sum-est trick log Σ; e:: max, Ijj + log(2 2Ex - #:). Since the matrix A is small we compute its inverse using double precision to improve stability. Surprisingly, we did not observe any benefit from using the Cholesky decomposition not calculating A-1 and instead solving the equation B = AX for X. We furthermore precompute W A- VV to avoid unnecessary operations.
We use layers and an embedding size Hs 32 for GTN. The MLPs use single hidden layer, biases and LeakyReLU non-linearities. The single-head MLP uses an output size of Hs, match = HN and a hidden embedding size of 4HN, i.e. the same as the concatenated node embedding, and the multi-head MLP uses hidden embedding size of HN. To stabilize initial training we scale the node embeddings by
directly before calculating OT. d denotes the average graph distance in the training set, ñ the average number of nodes per graph, and Hl,, match the matching embedding size, i.e. 32 for single-head and 128 for multi-head OT.
Rv mack
J GRAPH DATASET GENERATION AND EXPERIMENTAL DETAILS
The dataset statistics are summarized in Table 5
Each dataset contains the distances between all graph pairs in each split, i.c. 10 296 and 1 1128 distances for preferential attachment. The AIDS dataset was generated by randomly sampling graphs with at most 30 nodes from the original AIDS dataset (Riesen Bunke 2008). Since not all node types are present in the training set and our choice of GED is permutation innmuraaioonnnvariintt
invariant w.r.t. types, we permuted the node types so that there are no previously unseen types in the validation and test sets. For the preferential attachment datasets we first generated 12, 4, and 4 undirected "seed" graphs (for train, val, and test) via the initial attractiveness model with randomly chosen parameters: 1 to 5 initial nodes, initial attractiveness of 0 to and 1/2n and 3/2m total nodes, where is the average number of nodes (20, 200. 2000, and 20 000). We then randomly label every node (and edge) in these graphs uniformly. To obtain the remaining graphs we edit the "seed" graphs between n/40 and n ni/20 times by randomly adding, type editing, removing nodes
20

---

Under review as a conference paper at ICLR 2021
Table Hyperparameters for the Linux dataset.
Ir
batchsize
layers
emb. size
La reg.
Abase
SiamMPNN
1 x 10
256
3
32
x 10
GMN
1 x 10
20
3
64
GTN, head
0.01
1000
32
10
1.0
8 OT heads
0.01
1000
3
32
x 10
1.0
Balanced OT
0.01
1000
3
32
x 10
2.0
Table Hyperparameters for the AIDS dataset.
Ir
batchsize
layers
emb. size
L. reg.
Abase
SiamMPNN
X 10
256
32
5 x 10
SimGNN
X 10-3
32
0.01
GMN
10-2
128
32
0
GTN, head
0.01
100
3
32
5 10-3
0.1
8 OT heads
0.01
100
3
32
5 10-3
0.075
Balanced OT
0.01
100
3
32
5 10-3
0.1
Nyström
0.015
100
32
5 10-3
0.2
Multiscale
0,015
100
3
32
5 10-3
0.2
Sparse OT
0.015
100
3
32
5 10-3
0.2
LCN-OT
0.015
100
32
10-3
0.2
and edges. Editing nodes and edges 4x and adding/deleting edges as likely as adding/deleting nodes. Most of these numbers were chosen arbitrarily, aiming to achieve somewhat reasonable dataset and process. We found that the process of first generating seed graphs and subsequently editing these is crucial for obtaining meaningfully structured data to learn from. For the GED choose an edit cost of for changing a node or edge type and for adding or deleting a node or an edge.
We represent node and edge types as one-hot vectors. We train all models except SiamMPNN (which uses SGD) and GTN Linux with the Adam optimizer and mean squared error (MSE) loss for up to 300 epochs and reduce the learning rate by a factor of 10 every 100 steps. On Linux we train for up 1000 epochs and reduce the learning rate by a factor of 2 every 100 steps. We use the parameters from the best epoch based on the validation set. We choose hyperparameters for all models using multiple steps of grid search on the validation set, see Tables 6 to for the final values. We use the originally published result of SimGNN on Linux and thus don i provide its hyperparameters. GTN uses 500 Sinkhorn iterations. We obtain the final entropy regularization parameter from Anase via where d denotes the average graph distance and the average number nodes The factor d/ñ serves to estimate the embedding distance scale and
_ 1_ Abuse n. To per graph in the training // log n counteracts the entropy scaling with n log n. Note that the entropy regularization parameter was small, but always far from O, which shows that entropy regularization actually has positive
Table 8: Hyperparameters for the preferential attachment GED dataset.
Ir
batchsize
layers
emb. size
L2 reg.
Abuse
SiamMPNN
1 10
256
64
1 X 10
SimGNN
x 10
32
0
GMN
1 x 10
20
3
64
0
GTN, head
0.01
100
3
32
5 x 10-1
0.2
OT heads
0.01
100
32
5 10-3
0.075
Balanced OT
0.01
100
32
5 10-4
0.2
Nyström
0.02
100
3
32
5 10
0.2
Multiscale
0.02
100
3
32
5 10-5
0.2
Sparse OT
0.02
100
32
10-5
0.2
LCN-OT
0.02
100
3
32
5 10-5
0.2
21

---

Under review as a conference paper at ICLR 2021
Table 9: Runtimes (ms) of Sinkhorn approximations for EN-DE embeddings at different dataset sizes. Full Sinkhorn scales quadratically, while all approximationes scale at most linearly with the size. Sparse approximations are 2-4x faster than low-rank approximations, and factored OT multiple times slower due to its iterative refinement scheme. Note that similarity matrix computation time (K)
primarily depends on LSH/Nyström method, not OT approximation.
N
10000
N
20000
N
50000
K
OT
K
OT
K
OT
Full Sinkhorn
2950
760
OOM
OOM
Factored OT
29
809
32
1016
55
3673
Multiscale OT
90
48
193
61
521
126
Nyström Skh.
29
135
41
281
79
683
Sparse Skh.
42
46
84
68
220
137
ACN-Sinkhorn
101
116
242
205
642
624
10000
Multse. OT Nys. Skh.
300
Time per epoch (s)
Runtime (ms)
Sparse Skh.
LCN-Skh.
1000
200
100
100
Full
LCN
10
0
100
200
1000
10000
Neighbors + landmarks
Avg. graph size
Figure Runtime scales linearly with the number of neighbors/landmate for all relevant Sinkhorn approximation methods.
Figure Log-log runtime per epoch for GTN with full Sinkhorn and LCN-Sinkhorn. LCN- Sinkhorn scales almost linearly with graph size while sustaining similar accuracy.
effect on learning. On the pref. 200 dataset no La regularization, Amase 0.5, and a batch size of 200. For pref. att. 2k we use Abase and a batch size of 20 for full Sinkhorn and 100 for LCN-OT. For pref. 20k we use Abase 50 and batch size of 4. Abase scales with graph size due to normalization of the PM kernel.
For LCN-OT we roughly 10 neighbors for LSH (20 k-means clusters) and 10 k-means landmarks for Nyström on pref. att. 200. We double these numbers for pure Nyström Sinkhorn, sparse OT, and multiscale OT. For pref. att. 2k we use around 15 neighbors (10 20 hierarchical clusters) and 15 landmarks and for pref. att. 20k we use roughly 30 neighbors (10 10 10 hierarchical clusters) and 20 landmarks. The number of neighbors for the 20k dataset higher and strongly varies per iteration due to the unbalanced nature of hierarchical k-means. This increase in neighbors and landmarks and Py PyTorch's 's missing support for ragged tensors largely explains LCN-OT's deviation from perfectly linear runtime scaling.
We perform all runtime measurements on a compute node using one Nvidia GeForce GTX 1080 Ti, two Intel Xeon E5-2630 v4, and 256GB RAM.
K RUNTIMES
Table 9 compares the runtime of the full Sinkhorn distance with different approximation methods using 40 neighbors/landmarks. We separate the computation approximate K from the optimal transport computation (Sinkhorn iterations), since the former primarily depends on the LSH and Nyström methods we choose. We observe a 2-4x speed difference between sparse (multiscale OT and sparse Sinkhorn) and low-rank approximations (Nyström Sinkhorn and LCN-Sinkhorn), while factored OT multiple times slower due to its iterative refinement scheme. In Fig. 5 we observe that this runtime gap stays constant independent of the number of neighbors/landmarks i.e. the relative
22

---

Under review as a conference paper at ICLR 2021
Fact. OT
2.0
2.5
Multse. OT
Nys. Skh,
1/00
0.0
Spurse Skh.
LLN. Skhh-Skh.
Rel. emr.
1.0
0.5
-5.0
100
200
100
200
0.0
10
10
10-1
10
10
Runtime (ms)
Neighbors + landmarks
Figure 7: Sinkhorn distance ap- Figure 8: Sinkhorn distance ap- Figure 9: Sinkhorn distance proximation at different runtimes proximation for varying compu- approximation for varying en- (varied via the number of neigh- tational budget. The dashed tropy regularization À at constant bors/landmarks), The dashed line denotes the true Sinkhorn computational budget. Sparse line denotes the true Sinkhorn distance. Sparse Sinkhorn Sinkhorn performs best for low distance. Sparse Sinkhorn con- mostly performs best, with LCN- X, LCN-Sinkhorn for moderate sistently performs best. The Sinkhorn coming in second. Fac- and high À and factored OT for row indicates factored OT results tored OT performs well with very high A. far outside the depicted range.
few landmarks.
difference decreases as we increase the number of neighbors/landmarks. This gap could either be due to details in low-level CUDA implementations and hardware or the fact that low-rank approximations require 2x many multiplications for the same number of neighbors/landmarks. In either case, both Table 9 and Fig. show that the runtimes of all approximations scale linearly both in the dataset size and the number of neighbors and landmarks, while full Sinkhorn scales quadratically.
We furthermore investigate whether GTN with approximate Sinkhorn indeed scales log-linearly with the graph size by generating preferential attachment graphs with 200, 2000, and 20 000 nodes (+-50 %). We use the Pyramid matching (PM) kernel (Nikolentzos et al. 2017) as prediction target. Fig. 6 shows that the runtime of LCCN-Sinkhorn scales almost linearly (dashed line) and regular full Sinkhorn quadraticly (dash-dotted line) with the number of nodes, despite both achieving similar accuracy and LCN using slightly more neighbors and landmarks on larger graphs sustain good accuracy. Full Sinkhorn went out of memory for the largest graphs.
L DISTANCE APPROXIMATION
Figs. and show that for the chosen À = 0.05 sparse Sinkhorn offers the best trade-off between computational budget and distance approximation, with LCN LCN- Sinkhorn and multiscale OT coming in second. Factored OT is again multiple times slower than the other methods and thus not included in Fig. 7 Note that d can be negative due to the entropy offset. This picture changes as we increase the regularization. For higher regularizations CNN-Sinhoom is the most precise constant computational budget (number of neighbors/landnarrke as shown in Fig. Note that the crossover points this figure roughly coincide with those in Fig. 4 Keep in mind that in most cases the OT plan is more important than the raw distance approximation, since it determines the training gradient and tasks like embedding alignment don't 't use the distance at all. This becomes evident in the fact that sparse Sinkhorn achieves a better distance approximation than LCN-Sinkhorn but performs worse in both downstream tasks investigated in [Sec. 7
23


Reviews:
[{"review_id": "dbtvdqdH5cT", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Proposal:\n- First log-linear time algorithms for entropy-regularized OT that work for complex \n  real-world tasks using high-dimensional spaces with little to no loss in accuracy\n  ... many claims in one statement\n- Locally Corrected Nyström\n  ... this would deserve a single paper - just to proof everything is still fine and valid\n  in particular to alternatives\n- In the way this paper is written I am positive it gets accepted \n  (because it fits the writing style of the more recent papers in the field)\n  and contains sufficient novelty\n  ... but I always wonder if good (published) science originates from clarity (or confusion)\n\ncomments:\n- sorry but the text is very hard to follow !\n- 'Optimal transport is concerned with the problem' (on p.3) - I think some introductive\n  work may not harm in the first page\n- the reader is thrown up by terms and references - in my view more confusing than enlighting\n  --> it may not harm to add some brief explainations of terms (from Cuturi:)\n  'A transport plan is a flow on that graph satisfying source (a i flowing out of each node i) \n   and sink (b j flowing into each node j 0 ) ... --- which is simple an optimal flow\n  in a graph ... I am not sure why we not simple can call it like this but need to come up with\n  new terms\n- The paper is written (following the very strange title\n   ... although Cuturi did the same it would be nice if we can stop having marketing titles\n  but focus a bit on science again ... in particular in 10 years many things proposed nowadays\n  are not lightspeed or warpspeed anymore) like providing an all issues solving theory \n  --> this does not improve the readability of the paper\n  For example Eq.2 what is the 'meaning' of (s) and (t) -- I have an idea but it is not written there\n- it is hard to proofread and verify a paper if it is written with the objective to confuse the\n  reviewer ;-)\n- widely incremental work by combing some known ideas (Nystreom, LSH, ...) - with a lot of addon\n  theory which is probably correct but not very clear in the presentation\n- 'Since the Nyström method is a low-rank approximation it only accounts for the\nglobal structure of the kernel matrix K and not for the local structure around each point x.'\n  - this is actually wrong (!) - if the landmarks are indepentent and the number of landmarks\n  aligns with the rank of the matrix - Nystroem will provide a perfect (!) reconstruction\n\t--> there is a lot of work on the approximation bound of Nystroem (and related methods) - see\n  e.g. work by Dhillon\n- where is the definition of 'sparse approximation K^sp' used in Eq 3? \n- how precisely does \\bar{P} (after Eq 3) link to the part around Eq 2? - is the Kernel K_{ij} used\n  here as well and / or where is the actually input data Kernel matrix \n- in Eq 1 what should be a cost function here and how do you obtain C_{ij}?\n- Ok Eq 4 is an actual proposal by balancing (and joining) sinkhorn and Nystroem in one distance formulation\n  and it would not harm to motivate why and where you need such a distance in advance \n  (problem statement --> solutions --> particular strategy --> outlined proposal --> evaluation + proofs)\n- 'Most modern ML models are trained using backpropagation' - lets rephrase it as: nowadays neural network\n  approaches are trained by backprop ... there are many other methods which are not at all trained by backprob\n  for good reasons \n- '... Usually we want to learn embeddings which act as point sets X p and X q and therefore need gradients' \n  - well yes, if we stick on neural networks we need vectorial inputs and hence are often looking for (costly)\n  embeddings - if we do not use NN we may not have this problem (but others)\n- 'We can either estimate these gradients via automatic differentiation' - this is in general the more costly\n  way to do things and I am happy to see that explicit derivations are given\n- regarding table 1 --> before you come with numbers (where you measured something) it would be good to specify\n  details of your scenario (which are omitted before) - in particular which data, which cost function, which parameters \n  a.s.o. -- and although I understand that you like your method most it would still be good to provide some\n  oldfashion baseline (and not - not from sinkhorn)\n- 'We propose the graph transport network (GTN) to evaluate approximate Sinkhorn and enhanced optimal transport and advance the state of the art on this task.' --- fantastic on page 6 you actually outline a more userfriendly motivation\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "XGprfy8b9aw", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper studies how to approximate Sinkhorn computation using more efficient kernel matrix representations (low-rank approach + LSH based sparse approach). Neither of both ideas is completely new, the authors studies a combination of them that hasn't been explored in the literature, and use the proposed tech in three applications: ranking, embedding alignment, graph distance regression. \n\nPro:\n\n- Authors have attempted multiple applications to prove the effectiveness with quantitative metrics. The main contribution seems the empirical validations. \n\nCon:\n\n- I think the presentation is a bit out of focus. Some sections can be left for appendix (e.g., Backpropagation in Sec 4. and Sec 5.) since some are either very standard in the literature or not really solidly experimented. Since I consider the main contribution to be the empirical evidences, the space should left for more details on those empirical experiments (for reproducibility purpose).\n\n- the paper's idea is not particularly innovative. (yet it is not the sole reason for my scoring). \n\n- Some relevant works on low-rank ideas have not been compared/cited. e.g.  \n\nForrow, Aden, et al. \"Statistical optimal transport via factored couplings.\" The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.\n\n- The writing can be improved. A lot citations/references are not particularly relevant to what this paper is about and make some parts not enough self-explained. \n\n----\nThe authors addressed my concerns and I raised my evaluation ratings. \n\n\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "RMShAVkNfCr", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\nThe paper considers the problem of approximating Sinkhorn divergence and corresponding transportation plan by combining low-rank and sparse approximation for the Sinkhorn kernel and using Nystrom iterations as a substitute for Sinkhorn's iterations. The corresponding approach is amenable to differentiation and can be used as a building block in different architectures. Numerical experiments in several settings are performed to compare  the proposed approach with existing ones and demonstrate its scalability.\n\nEvaluation:\nI believe the proposed framework is a valuable contribution in terms of practical performance and wide list of applications where OT could not be used before because of the high computational cost. So, I would recommend accepting this paper. \n\nPros:\n1. High scalability of the proposed approach and linear up to log factors in dimension complexity.\n2. Flexibility of the framework due to a combination of sparse and low-rank approximations, which are complementary to each other.\n\nCons:\n1. Some parts of the paper seem to be not clear for a general audience.\n\na. First page. $n$ is undefined.\n\nb. First paragraph of Sect. 2. What is \"set of embeddings\"?\n\nc. Last but one paragraph on p.2. $d$ is not defined.\n\nd. In (1) $F$ stands for the Frobenius product, does it?\n\ne. Proposition 1. $N$ is not defined.\n\nf. First paragraph of Sect. 5. What is \"OT with multiple heads\"?\n\ng. What is meant as embedding?\n\nh. In the experiments, what is used as the cost function to define the Sinkhorn kernel $K$? If this is an $L_2$ distance, the convolutions can be used to accelerate the standard Sinkhorn and it would be nice to see the comparisons with convolutional Sinkhorn, which is also log linear.\n\ni. Appendix A. $B,r,b$ are not defined when they are first used.\n\nj. In (17), how was the last equality obtained?\n\nk. In (19), (20), how were the first equalities obtained?\n\nl. Appendix E, first paragraph. What is \"log-sum-exp trick\"?\n\nm. Appendix G. What is \"similarity matrix\"?\n\n2. As far as I understood, the proposed approach is not amenable to parallel computations on GPU as opposed to standard Sinkhorn.\n\n\nMinor comments\n1. Maybe it is too strong to state in the abstract that this is the first log-linear time algorithm given that when the Sinkhorn kernel corresponds to a convolution, the Sinkhorn's algorithm is log-linear by using the FFT.\n2. Bibliographical note. (Altschuler et al., 2017) did not show $1/\\varepsilon^2$ bound for the Greenkhorn. Their bound for Greenkhorn is the same $1/\\varepsilon^3$ as for the Sinkhorn. The bound for Sinkhorn was improved to $1/\\varepsilon^2$ in http://proceedings.mlr.press/v80/dvurechensky18a.html and the bound for Greenkhorn was improved to $1/\\varepsilon^2$ in http://proceedings.mlr.press/v97/lin19a.html.\n3. Bibliographical note. Quadratic regularization for OT was proposed in https://arxiv.org/abs/1704.08200.\n4. Appendix A. I believe that in this framework a general value of the regularization parameter $\\lambda$ is used. If it is the case, then the number of Sinkhorn iterations to find an $\\varepsilon$-solution to the regularized problem is $1/(\\varepsilon \\lambda)$. This follows from http://proceedings.mlr.press/v80/dvurechensky18a.html Theorem 1 and an estimate for $R$ in Lemma 1. The bound $1/\\varepsilon^2$ corresponds to finding an $\\varepsilon$-solution for the non-regularized problem. In this case one has to set $\\lambda=\\varepsilon/(4 \\ln n)$, which may be too small.", "extra_scores": {}, "extra_text": {}}, {"review_id": "7E8KZoQ0t-C", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Overall, I found this paper interesting, and I think it does address a relevant problem for the community.\n\n\nI have been playing with Nystrom approximations myself and I know the results are a bit disappointing, but it is grounded on strong theory. Then, it is pretty much welcome that attempts to patch the problems of Nystroms are made, so Optimal Transport becomes more scalable.\n\nThe reason for why my judgement is below acceptance is that I believe both theory and results altogether are not strong enough so they live up to what is promised in the abstract. Typically when new methods are proposed with the promise of bridging a gap and solving a relevant problem, I hope they will have a thorough theoretical justification, or the results will be compelling. None of this convincing enough here.\n\n1)On the theoretical side, I missed a convergence analysis, as the one in the Nystrom method. The theory side focuses on some derivative calculations, but I would love to see how the interplay between sparsity, Nystrom method and LSH lead to better convergence. I acknowledge this can be hard to do, though. Without having the theory it is hard to understand whether any reported experimental result is a consequence of choosing particularly good example for their method.\n\n2)I found the experimental/methodological side was a bit disconnected from the rest; the paper contains  several vignettes about some applications/improvements in the practical side, but when reading that felt like belonging to other paper. I recommend authors work in creating a more coherent story\n2a) I found the discussion on Multi-Head OT unjustified and even a bit misleading. The Authors refer to some NLP papers, like arguing OT plays a role there. But none of these papers have any OT at all. The analogy of  \"softmax for rows\" is not convincing as this is simply a softmax applied many times. There is a world of difference between that and the result of sinkhorn algorithm, but the narrative seems to downplay the actual difference between them. I recommend the authors elaborate more on this connection, because otherwise it is hard to follow (and the subsequent results).\n2b)Results on translation seem impressive (Table 2), but raise a concern. Why would your method outperform Sinkhorn if it is only approximation? is it perhaps the result of randomness? since an explanation of this phenomenon is missing I am led to believing. Authors should improve the exposition of the baseline \"original\". Why does full Sinkhorn does better than original?. In summary, I think authors should improve the discussion about the validity/significance of their empirical results, highlighting the regimes when they are supposed to express and when they are not.\n2c)The main figure is Fig 2. I recommend authors build on this and expand those results so it is clear when their method is better and when it is not", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors propose to approximate the kernel matrix used in the Sinkhorn algorithm by a combination of sparse + low rank approximation. To do so, the authors propose to compute a low rank approximation of a sparsified (thresholded below a certain value to be 0) kernel matrix using Nyström, and then correct it by adding back the true entries at non-sparse entries, after removing those obtained from the approximation. This results in a matrix whose application then results in sparse + low-rank.

The first version of the paper contained mostly experimental evidence, which was deemed a bit short by some reviewers.The authors have added theoretical material on the way. Although I believe these are worthy additions, as AC, I do not feel comfortable accepting the paper as of now, because I believe these additions were not properly reviewed. I understand this must be disappointing for the authors, who have sprinted to add new content during the rebuttal phase, but I hope they agree that the rebuttal process is not here to handle entirely new sections, but rather to improve existing parts. In particular, that section should be reviewed by authors knowledgeable on low rank kernel factorization, something I did not see in the pool of reviewers. I also believe the paper still has a few shortcomings. Taken together, I therefore recommend a re-submission.

ideas to improve the paper

- the authors claim to use Nyström on a sparsified matrix (see eq. 4). The sparsified kernel is no longer positive definite. I would like the authors to comment on this. I understand Nyström could be used naively without any psd-ness guarantees, but I think a heads-up is needed.There are, furthermore, several local/global factorizations of kernel matrices available out there (e.g. MEKA, https://www.jmlr.org/papers/v18/15-025.html), the main difference here being that the product by such approximation must be guaranteed to be positive for it to work in the Sinkhorn algorithm. I would expect that bounds in expectation to break down sometimes, and therefore result in "catastrophic" failures (i.e. nan's). I think that an algorithm that claims to improve or replace another one, and which has such blind spots, needs such additional experiments (I have read the Limitations section in Appendix B, something more precise would improve the paper). I understand these were not part of the original Nyström paper for Sinkhorn, but since this is an increment over that previous work, therefore lacking a bit its originality, more knowledge needs to be contributed.

- For instance, since the authors write an entire paragraph on this (Appendix B), I am surprised that there is not direct mention to the fact that a sparse sinkhorn may simply *not* converge, because it may not satisfy the fully indecomposable property required of matrices for Sinkhorn's algorithm to converge. 

- i dont think that users have the various identities (14,15) in mind when they think about "backpropagating" through Sinkhorn. What is typically needed is to compute the differentiable properties of the regularized OT matric and/or of the regularized OT cost w.r.t. *point locations* (i.e. x_i). The statement "LCN-Sinkhorn works flawlessly with automatic backpropagation" is misleading in the sense that it ignores that problem altogether. Since so many extensions of OT today relay on that differentiability, the section, as it is written now, is problematic.

- several methods claim to be faster of more efficient than Sinkhorn to solve OT. Either these methods display faster theoretical convergence (e.g. by using acceleration) or display faster practical convergence (e.g. heavy ball variants) using synthetic, controlled datasets. Using synthetic data helps exhibit highlight relevant regimes for regularization parameters, including those where LSE Sinkhorn may converge but LCN does not work, or vice-versa. I understand that the authors' wanted to use real data, but it would be great to clarify whether that setup was used because LCN works better there (in which case this becomes more of a paper at the intersection of OT and word embeddings) or because this happened to be the first and only example the authors thought of.
Decision: Reject