Year: 2020

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.

heading: 1 INTRODUCTION
text: Measuring the distance between two distributions or sets of objects is a 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 a source distribution (e.g. a pile of earth) to a 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., 2015), set matching (Wang et 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; Papadakis et al., 2014; 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 high-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-regularized optimal transport: Sparse Sinkhorn and LCNSinkhorn, 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-regularized 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 LCN, a general similarity matrix approximation that fuses local (sparse) and global (low-rank) approximations, allowing us to simultaneously capture both kinds of behavior. LCN-Sinkhorn thus fuses sparse Sinkhorn and Nyström-Sinkhorn (Altschuler et al., 2019). Both sparse Sinkhorn and LCN-Sinkhorn run in log-linear 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 a wide range of regularization parameters and computational
budgets (as e.g. demonstrated in Fig. 1). 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 et 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, we develop the graph transport network (GTN), which combines graph neural networks (GNNs) with optimal transport, and further improve it 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:
• Locally Corrected Nyström (LCN), a flexible, log-linear time approximation for similarity matrices, leveraging both local (sparse) and global (low-rank) approximations.
• Entropy-regularized optimal transport (a.k.a. Sinkhorn distance) with log-linear runtime via sparse Sinkhorn and LCN-Sinkhorn. These are the first log-linear approximations that are stable enough to substitute full entropy-regularized OT in models that leverage high-dimensional spaces.
• The graph transport network (GTN), which combines a graph neural network (GNN) with multihead unbalanced LCN-Sinkhorn. GTN both sets the state of the art on graph distance regression and still scales log-linearly in the number of nodes.

heading: 2 SPARSE SINKHORN
text: Entropy-regularized 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 ∈ Rn and q ∈ Rm supported on two sets of points Xp = {xp1, . . . ,xpn} and Xq = {xq1, . . . ,xqm} in Rd and the cost function c : Rd × Rd → R (e.g. the squared L2 distance) giving rise to the cost matrix Cij = c(xpi,xqi) we aim to find the Sinkhorn distance dλc and the associated optimal transport plan P̄ (Cuturi, 2013)
dλc = min P 〈P ,C〉F − λH(P ), s.t. P1m = p, P T1n = q, (1)
with the Frobenius inner product 〈., .〉F and the entropy H(P ) = − ∑n i=1 ∑m j=1Pij logPij . Note that dλc includes the entropy and can thus be negative, while Cuturi (2013) originally used d 1/λ Cuturi,c = 〈P̄ ,C〉F. This optimization problem can be solved by finding the vectors s̄ and t̄ that normalize the columns and rows of the matrix P̄ = diag(s̄)K diag(t̄) with the similarity matrixKij = e− Cij λ , so that P̄1m = p and P̄ T1n = q. This is usually achieved via the Sinkhorn algorithm, which initializes the normalization vectors as s(1) = 1n and t(1) = 1m and then updates them alternatingly via
s(i) = p (Kt(i−1)), t(i) = q (KTs(i)) (2) until convergence, where denotes elementwise division. Sparse Sinkhorn. The Sinkhorn algorithm is faster than non-regularized EMD algorithms, which run in O(n2m log n log(nmax(C))) (Tarjan, 1997). However, its computational cost is still quadratic
in time, i.e. O(nm), which is prohibitively expensive for large n and m. We propose to overcome this by observing that the matrix K, and hence also P̄ , is negligibly small everywhere except at each point’s closest neighbors because of the exponential used inK’s computation. We propose to leverage this by approximating C via the sparse matrix Csp, where
Cspij = { Cij if xpi and xqj are “near”, ∞ otherwise. (3)
Ksp and P̄ sp follow according to the definitions of K and P̄ . In this work we primarily consider neighbors with distance lower than r1 as “near”. Finding such neighbors can be efficiently solved via locality sensitive hashing (LSH) on Xp ∪Xq. 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 r1 are put into the same bucket with probability at least p1, while those beyond some distance r2 = c · r1 with c > 1 are put into the same bucket with probability at most p2 p1. 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. O(n log n) for n ≈ m (see App. A and App. K for details). Unfortunately, LSH can fail when e.g. the cost between pairs is very similar (see App. B). However, we can alleviate these limitations by fusingKsp with the Nyström approximation.

heading: 3 LOCALLY CORRECTED NYSTRÖM AND LCN-SINKHORN
text: Nyström method. The Nyström method is a popular way of approximating similarity matrices that provides performance guarantees for many important tasks (Williams & Seeger, 2001; Musco & Musco, 2017). It approximates a positive semi-definite (PSD) similarity matrixK via its low-rank decomposition KNys = UA−1V . Since the optimal decomposition via SVD is too expensive to compute, Nyström instead chooses a set of l landmarks L = {xl1, . . . ,xll} and obtains the matrices via Uij = k(xpi,xlj), Aij = k(xli,xlj), and Vij = k(xli,xqj), where k(x1,x2) is an arbitrary PSD kernel, e.g. k(x1,x2) = e− c(x1,x2)
λ for Sinkhorn. Common methods of choosing landmarks from Xp ∪Xq are uniform and ridge leverage score (RLS) sampling. We instead focus on k-means Nyström and sampling via k-means++, which we found to be significantly faster than recursive RLS sampling (Zhang et al., 2008) and perform better than both uniform and RLS sampling (see App. H).
Sparse vs. Nyström. Exponential kernels like the one used forK (e.g. the Gaussian kernel) typically have a reproducing kernel Hilbert space that is infinitely dimensional. The resulting Gram matrixK 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 x. As such, it is ill-suited for any moderately low entropy regularization parameter, where the transport matrix P̄ resembles a 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 as such not competing but rather complementary approaches.
Locally corrected Nyström. Since we know that the entries in our sparse approximation are exact, fusing this matrix with the Nyström method is rather straightforward. For all non-zero values in the sparse approximationKsp we first calculate the corresponding Nyström approximations, obtaining the sparse matrixKspNys. To obtain the locally corrected Nyström (LCN) approximation we remove these entries fromKNys and replace them with their exact values, i.e.
KLCN = KNys +K sp ∆ = KNys −K sp Nys +K sp. (4)
LCN-Sinkhorn. To obtain the approximate transport plan P̄LCN 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) viaKLCNt = U(A−1V t) +K sp ∆t, similarly to Altschuler et al. (2019). As a result we obtain the decomposition of P̄LCN = P̄Nys + P̄ sp ∆ = P̄U P̄W + P̄ sp − P̄ spNys and the approximate distance (using Lemma A from Altschuler et al. (2019))
dλLCN,c = λ ( sT P̄U P̄W1m + 1 T n P̄U P̄W t+ s T P̄ sp∆ 1m + 1 T n P̄ sp ∆ t ) . (5)
This approximation scales log-linearly with dataset size (see App. A and App. K for details). It allows us to smoothly move from Nyström-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.

heading: 4 THEORETICAL ANALYSIS
text: 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 ofK 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. C to G. Theorem 1. Let Xp and Xq have n samples that are uniformly distributed in a d-dimensional closed, locally Euclidean manifold with unit volume. Let furthermore Cij = ‖xpi − xqj‖2 and Kij = e
−Cij/λ. Let the l landmarks L be arranged optimally and regularly so that the expected L2 distance to the closest landmark is minimized. Denote R = 12 minx,y∈L,x6=y ‖x− y‖2. Assume that the sparse correctionKspij = Kij if and only if xqj is one of the k − 1 nearest neighbors of xpi, and that the distance to xpi’s k-nearest neighbor δk R. Then the expected maximum error in row i of the LCN approximationKLCN is
E[‖Ki,: −KLCN,i,:‖∞] = E[e−δk/λ]− E[KLCN,i,j ], (6) with j denoting the index of xpi’s k-nearest neighbor. Using the upper incomplete Gamma function Γ(., .) we can furthermore bound the second term by
e− √ dR/λ ≤ E[KLCN,i,j ] ≤ 2d(Γ(d)− Γ(d, 2R/λ)) (2R/λ)d(1 + e−2R/λ) +O(e−2 √ 3R/λ). (7)
The error in Eq. (6) is dominated by the first term since δk R. Note that R only decreases slowly with the number of landmarks since R ≥ ( (d/2)!l ) 1/d 1 2 √ π
(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 we obtain a 68 % improvement in the first term (d=32, λ=0.05, n=1000). This is even more pronounced in clustered data. Theorem 2. Let Xp, Xq ⊆ Rd be distributed inside the same c clusters with cluster centers xc. Let r be the maximum L2 distance of a point to its cluster center and D the minimum distance between two points from different clusters, with r D. Let each LSH bucket used for the sparse approximation Ksp cover at least one cluster. Let KNys use 1 ≤ l ≤ d and KLCN use l = 1 optimally distributed landmarks per cluster. Then the maximum error is
max ‖K −KNys‖∞ = 1− max ∆∈[0,r]
le−2 √ r2+ l−12l ∆ 2/λ 1 + (l − 1)e−∆/λ −O(e−D/λ), (8)
max ‖K −Ksp‖∞ = e−D/λ, (9) max ‖K −KLCN‖∞ = e−D/λ(1− e−2r/λ(2− e−2r/λ) +O(e−D/λ)). (10)
Since we can lower bound Eq. (8) by 1 − le−2r/λ − O(e−D/λ) we can conclude that the error in KNys is close to 1 for any reasonably large rλ (which is the maximum error possible). The errors in Ksp andKLCN on the other hand are vanishingly small, since r D. Moreover, these maximum approximation error improvements directly translate to 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 et al. (2019). Theorem 3 (Altschuler et al. (2019)). Let Xp, Xq ⊆ Rd have n samples. Denote ρ as the maximum distance between two samples. Let K̃ be an approximation of the similarity matrixK withKij = e−‖xpi−xqj‖2/λ and ‖K̃ −K‖∞ ≤ ε ′ 2 e −ρ/λ, where ε′ = min(1, ε
50(ρ+λ log λnε ) ). When performing
the Sinkhorn algorithm until ‖P̃1N − p‖1 + ‖P̃ T1N − q‖1 ≤ ε′/2, the resulting approximate transport plan P̃ and distance d̃λc are bounded by
|dλc − d̃λc̃ | ≤ ε, DKL(P̄ ‖P̃ ) ≤ ε/λ. (11)
Convergence rate. We next show that approximate Sinkhorn converges as fast as regular Sinkhorn by slightly adapting the convergence bound by Dvurechensky et al. (2018) to account for sparsity.
Theorem 4 (Dvurechensky et al. (2018)). Given the matrix K̃ ∈ Rn×n and p, q the Sinkhorn algorithm gives a transport plan satisfying ‖P̃1N − p‖1 + ‖P̃ T1N − q‖1 ≤ ε in iterations
k ≤ 2 + −4 ln(mini,j{K̃ij |K̃ij > 0}mini,j{pi, qj}) ε . (12)
Backpropagation. Efficient gradient computation is almost as important for modern deep learning models as the algorithm itself. These models usually aim at learning the embeddings in Xp and Xq and therefore need gradients w.r.t. the cost matrix C. We can estimate these either via automatic differentiation of the unrolled Sinkhorn iterations or 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 LSH neighbors, for which we use a simple straight-through estimator (Bengio et al., 2013)). 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. Proposition 1. The derivatives of the distances dλc and dλLCN,c (Eqs. (1) and (5)) and the optimal transport plan P̄ ∈ Rn×m w.r.t. the (decomposed) cost matrix C ∈ Rn×m in entropy-regularized OT and LCN-Sinkhorn are
∂dλc ∂C = P̄ , ∂P̄ij ∂Ckl = − 1 λ P̄ijδikδjl, (13)
∂dλLCN,c ∂U = −λs̄(Wt̄)T , ∂dλLCN,c ∂W = −λ(s̄TU)T t̄T , ∂dλLCN,c ∂ logKsp = −λP̄ sp, ∂dλLCN,c ∂ logKspNys = −λP̄ spNys,
(14)
∂P̄U,ij ∂Ukl = δikδjlsi, ∂P̄W,ij ∂Ukl = P̄ †U,ikskP̄W,lj , ∂P̄U,ij ∂Wkl = P̄U,iktlP̄ † W,lj ,
∂P̄W,ij ∂Wkl
= δikδjltj , ∂P̄ spij
∂ logKspkl = P̄ spij δikδjl,
∂P̄ spNys,ij ∂ logKspNys,kl = P̄ spNys,ijδikδjl,
(15)
with δij denoting the Kronecker delta and † the Moore-Penrose pseudoinverse. Using these decompositions we can backpropagate through LCN-Sinkhorn in time O((n+m)l2 + l3).

heading: 5 GRAPH TRANSPORT NETWORK
text: Graph distance learning. The ability to predict similarities or distances between graph-structured objects is 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., 2015), 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 of the art on this task.
Graph transport network. GTN first uses a Siamese graph neural network (GNN) to embed two graphs independently as sets of node embeddings. These embeddings are then matched using enhanced entropy-regularized optimal transport. Given an undirected graph G = (V, E), with node set V and edge set E , node attributes xi ∈ RHx and (optional) edge attributes ei,j ∈ RHe , with i, j ∈ V , we update the node embeddings in each GNN layer via
h (l) self,i = σ(W (l) nodeh (l−1) i + b (l)), (16)
h (l) i = h (l) self,i + ∑ j∈Ni η (l) i,jh (l) self,jWedgeei,j , (17)
with Ni denoting the neighborhood of node i, h(0)i = xi, h (l) i ∈ RHN for l ≥ 1, the bilinear layer Wedge ∈ RHN×HN×He , and the degree normalization η(1)i,j = 1 and η (l) i,j = 1/ √ degi degj for l > 1.
This choice of ηi,j allows our model to handle highly skewed degree distributions while still being able to represent node degrees. We found the choice of non-linearity σ not to be critical and chose a LeakyReLU. We do not use the bilinear layer Wedgeei,j if there are no edge attributes. We aggregate each layer’s node embeddings to obtain the final embedding of node i
hfinali = [h (1) self,i ‖h (1) i ‖h (2) i ‖ . . . ‖h (L) i ]. (18)
Having obtained the embedding sets Hfinal1 and H final 2 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 Hfinal1 and H final 2 . We obtain the embeddings for matching via h (0) i = MLP(h final i ) and obtain the final prediction via d = dλcwout + bout, with learnable wout and bout. All weights in GTN are trained end-to-end via backpropagation. For small graphs we use the full Sinkhorn distance and scale to large graphs by leveraging LCN-Sinkhorn. 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.
Learnable unbalanced OT. Since GTN regularly encounters graphs with disagreeing numbers of nodes it needs to be able to handle cases where ‖p‖1 6= ‖q‖1 or where not all nodes in one graph have a corresponding node in the other and thus P1m < p or P T1n < 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 et 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 C(p,ε)
C(ε,q) C(ε,ε)
] , C
(p,ε) ij =
{ ci,ε i = j
∞ i 6= j , C (ε,q) ij =
{ cε,j i = j
∞ i 6= j , C (ε,ε) ij = 0,
(19)
and adaptively computing the costs ci,ε, cε,j and cε,ε based on the input sets Xp and Xq. Using the BP matrix adds minor computational overhead since we only need to save the diagonals cp,ε and cε,q ofCp,ε andCε,q. We can then include the additional parts ofCBP in the Sinkhorn algorithm (Eq. (2)) via
KBPt = [ Kt̂+ cp,ε ť cε,q t̂+ 1Tn ť ] , KTBPs = [ KT ŝ+ cε,q š cp,ε ŝ+ 1Tmš ] , (20)
where t̂ denotes the upper and ť the lower part of the vector t. To calculate dλc we can decompose the transport plan PBP in the same way as CBP, with a single scalar for Pε,ε. For GTN we obtain the deletion cost via ci,ε = ‖α xpi‖2, with a learnable vector α ∈ Rd. Multi-head OT. Inspired by attention models (Vaswani et al., 2017) we further improve GTN by using multiple OT heads. UsingK heads means that we calculate OT in parallel forK separate sets of embeddings representing the same pair of objects and obtain a set of distances dλc ∈ RK . We can then transform these distances to a final distance prediction using a set of linear layers h(k)i = W
(k)hfinali for head k and obtain the final prediction via d = MLP(dλc ). Note that both learnable unbalanced OT and multi-head OT might be of independent interest.

heading: 6 RELATED WORK
text: 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 (Papadakis et al., 2014; Solomon et al., 2014), convolutions (Solomon et al., 2015), or embedding/hashing schemes (Indyk & Thaper, 2003; Andoni et al., 2008). In higher dimensions we can use tree-based algorithms (Backurs et al., 2020) or hashing schemes (Charikar, 2002), which are however limited to a previously fixed set of points Xp, Xq, on which only the distributions p and q change. For sets that change dynamically (e.g. during training) one common method of achieving log-linear runtime is a multiscale approximation of entropy-regularized OT (Schmitzer, 2019; Gerber & Maggioni, 2017). Tenetov et al. (2018) recently proposed using a low-rank approximation of the Sinkhorn similarity
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 Sec. 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 L2 distance as a cost function and are either unstable in convergence or prohibitively expensive for high-dimensional problems (O(nd3)) (Meng et al., 2019). Fast Sinkhorn. Another line of work has been pursuing accelerating entropy-regularized OT without changing its computational complexity w.r.t. the number of points. Original Sinkhorn requires O(1/ε2) iterations (Dvurechensky et al., 2018) and Jambulapati et al. (2019) recently proposed an algorithm that reduces them toO(1/ε). 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 an overview of cross-lingual word embedding models see Ruder et al. (2019). Unsupervised word embedding alignment was proposed by Conneau et al. (2018), with subsequent advances by Alvarez-Melis & Jaakkola (2018); 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 a single fixed-dimensional graph embedding (Bai et al., 2019; Li et al., 2019), or only use attention or some other strongly simplified variant of optimal transport (Bai et al., 2019; Riba et al., 2018; Li et al., 2019). Others break permutation invariance and are thus ill-suited for this task (Ktena et al., 2017; Bai et al., 2018). So far only approaches using a single graph embedding allow faster than quadratic scaling in the number of nodes. Compared to the Sinkhorn-based image model concurrently proposed by Wang et al. (2019) GTN uses no CNN or cross-graph attention, but an enhanced GNN and embedding aggregation scheme. OT has recently been proposed for graph kernels (Maretic et al., 2019; Vayer et al., 2019), which can (to a limited extent) be used for graph matching, but not for distance learning.

heading: 7 EXPERIMENTS
text: Approximating Sinkhorn. We start by directly investigating different Sinkhorn approximations. To do so we compute entropy-regularized OT on pairs of 10 000 word embeddings from Conneau et al. (2018), which we preprocess with Wasserstein Procrustes alignment in order to 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 λ = 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-Sinkhorn and factored OT (Forrow et al., 2019) obtain distances that are significantly closer to the true dλc than Multiscale OT and Nyström-Sinkhorn. Furthermore, the transport plan computed by sparse Sinkhorn and LCN-Sinkhorn show both a PCC and IoU that are around twice as high as Multiscale OT, while Nyström-Sinkhorn and factored OT exhibit almost no correlation. LCN-Sinkhorn performs especially well in this regard. This is also
evident in Fig. 1, which shows how the 104 × 104 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 quality. Factored OT exhibits a runtime 2 to 10 times longer than the competition due to its iterative refinement scheme. LCN-Sinkhorn performs best for use cases with constrained memory (few neighbors/landmarks), as shown in Fig. 3. The number of neighbors and landmarks directly determines memory usage and is linearly proportional to the runtime (see App. K). Fig. 9 shows that sparse Sinkhorn performs best for low regularizations, where LCN-Sinkhorn fails due to 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 be calculated). Interestingly, all approximations except factored OT seem to fail at high λ. We defer analogously discussing the distance approximation to App. L. All approximations scale linearly both in the number of neighbors/landmarks and dataset size, as 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 matrixR ∈ Rd×d 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 Conneau et al. (2018) by migrating Grave et al. (2019)’s implementation to PyTorch. The only change we make is using the full set of 20 000 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 use as many landmarks and neighbors as can fit into GPU memory and finetune both methods. Table 2 shows that using full Sinkhorn yields a significant improvement in accuracy on this task
Table 3: RMSE for GED regression across 3 runs and the targets’ standard deviation σ. GTN outperforms previous models by 48 %.
Linux AIDS30 Pref. att. σ 0.184 16.2 48.3 SiamMPNN 0.090(7) 13.8(3) 12.1(6) SimGNN 0.039 4.5(3) 8.3(1.4) GMN 0.015() 10.3(6) 7.8(3) GTN, 1 head 0.022(1) 3.7(1) 4.5(3) 8 OT heads 0.012(1) 3.2(1) 3.6(2) Balanced OT 0.034(1) 15.3(1) 27.4(9)
Table 4: RMSE for graph distance regression across 3 runs. Using LCN-Sinkhorn with GTN increases the error by only 10 % and allows log-linear scaling.
GED PM [10−2] AIDS30 Pref. att. Pref. att. 200
σ 16.2 48.3 10.2 Full Sinkhorn 3.7(1) 4.5(3) 1.27(6) Nyström Skh. 3.6(3) 6.2(6) 2.43(7) Multiscale OT 11.2(3) 27.4(5.4) 6.71(44) Sparse Skh. 44.0(30.4) 40.7(8.1) 7.57(1.09) LCN-Skh. 4.0(1) 5.1(4) 1.41(15)
compared to the original approach of performing Sinkhorn on randomly sampled subsets of embeddings (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 et 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 2 new datasets by computing the exact GED using the method by Lerouge et al. (2017) on small graphs (≤ 30 nodes) from the AIDS dataset (Riesen & Bunke, 2008) and a 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 on the validation set via a grid search. For more details see App. H to J. 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 we next ask whether we can sustain its performance 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 4 shows that both sparse Sinkhorn and the multiscale method using 4 (expected) neighbors 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 an overall low-rank structure with Nyström acting as regularization. However, it does not perform as well on the other two datasets. Using LCN-Sinkhorn with 2 neighbors and landmarks works well on all three datasets, with an RMSE increased by only 10 % compared to full GTN. App. K furthermore shows that GTN with LCN-Sinkhorn indeed scales linearly in the number of nodes across multiple orders of magnitude. This model thus allows to perform graph matching and distance learning on graphs that are considered large even for simple node-level tasks (20 000 nodes).

heading: 8 CONCLUSION
text: Locality sensitive hashing (LSH) and the novel locally corrected Nyström (LCN) method enable fast and accurate approximations of entropy-regularized OT with log-linear runtime: Sparse Sinkhorn and LCN-Sinkhorn. The graph transport network (GTN) is one example for such a 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-linearly 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.

heading: A COMPLEXITY ANALYSIS
text: Sparse Sinkhorn. A common way of achieving a high p1 and low p2 in LSH is via the AND-OR construction. In this scheme we calculate B · r hash functions, divided into B sets (hash bands) of r 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 asO((n+m)dBbr) for the hash functions we consider. As expected, for the tasks and hash functions we investigated we obtain approximately m/br and n/br neighbors, with br hash buckets per band. Using this we can fix the number of neighbors to a small, constant β in expectation with br = min(n,m)/β. We thus obtain a sparse cost matrix Csp with O(max(n,m)β) non-infinite values and can calculate s and t in linear time O(Nsink max(n,m)β), where Nsink ≤ 2 + −4 ln(mini,j{K̃ij |K̃ij>0}mini,j{pi,qj})ε (see Theorem 4) denotes the number of Sinkhorn iterations. Calculating the hash buckets with r = log min(n,m)−log βlog b takes O((n+m)dBb(log min(n,m)− log β)/ 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 n) for n ≈ m. LCN-Sinkhorn. Both choosing landmarks via k-means++ sampling and via k-means with a fixed number of iterations have the same runtime complexity of O((n+m)ld). PrecomputingW can be done in time O(nl2 + l3). The low-rank part of updating the vectors s and t can be computed in O(nl + l2 + lm), with l chosen constant, i.e. independently of n and m. Since sparse Sinkhorn with LSH has a log-linear runtime we again obtain log-linear overall runtime for LCN-Sinkhorn.

heading: B LIMITATIONS
text: 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 a sparse matrix. However, it can fail if the cost between pairs is very similar or the regularization is very high, if the dataset contains many hubs, i.e. points with a large number of neighbors, or if the distributions p or q are spread very unevenly. Furthermore, sparse Sinkhorn can be too unstable to train a model from scratch, since randomly initialized embeddings often have no close neighbors (see Sec. 7). LCN-Sinkhorn largely alleviates these limitations.
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 LCN-Sinkhorn to break 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 inverseA−1, 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 Ksp∆ . 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.

heading: C PROOF OF THEOREM 1
text: We first prove a lemma that will be useful later on. Lemma A. Let K̃ be the Nyström approximation of the similarity matrixKij = e−‖xi−xj‖2/λ. Let xi and xj be data points with equal L2 distance ri and rj to all l landmarks, which have the same distance ∆ to each other. Then
K̃ij = le−(ri+rj)/λ
1 + (l − 1)e−∆/λ (21)
Proof. The inter-landmark distance matrix is A = e−∆/λ1l×l + (1− e−∆/λ)Il, (22)
where 1l×l denotes the constant 1 matrix. Using the identity
(b1n×n + (a− b)In)−1 = −b
(a− b)(a+ (n− 1)b) 1n×n +
1
a− b In (23)
we compute
K̃ij = Ui,:A −1V:,j
= ( e−ri/λ e−ri/λ · · · )( −e−∆/λ (1− e−∆/λ)(1 + (l − 1)e−∆/λ) 1l×l + 1 1− e−∆/λ Il )e −rj/λ e−rj/λ
...  = e−(ri+rj)/λ
1− e−∆/λ
( −l2e−∆/λ
1 + (l − 1)e−∆/λ + l
) = e−(ri+rj)/λ
1− e−∆/λ l − le−∆/λ 1 + (l − 1)e−∆/λ
= le−(ri+rj)/λ
1 + (l − 1)e−∆/λ (24)
Now consider the error ‖Ki,: −KLCN,i,:‖∞. The k − 1 nearest neighbors are covered by the sparse correction and therefore the next nearest neighbor has distance δk. 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] ≥ EV (R)[r] = dd+1R. Because furthermore δk R, the error is dominated by the first term and the maximum error in row i is given by the k-nearest neighbor of i, denoted by j. Thus
E[‖Ki,: −KLCN,i,:‖∞] = E[Ki,j −KLCN,i,j ] = E[Ki,j ]− E[KLCN,i,j ] = E[e−δk/λ]− E[KLCN,i,j ] (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 Percus & Martin (1998))
E[e−δk/λ] = n!
(n− k)!(k − 1)!
∫ ((d/2)!)1/d√ π
0
e−r/λV (r)k−1(1− V (r))n−k dV (r) dr dr, (26)
with the volume of the d-ball
V (r) = πd/2rd
(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 r by √ dR/2, since this would be the furthest distance to the closest point in a d-dimensional grid. Any optimal arrangement minimizing E[miny∈L ‖x− y‖2 |x ∈ Xp] would be at least as good as a grid and thus have furthest distances as small or smaller than those in a grid. Thus,
E[KLCN,i,j ] (1) ≥ e−2r/λ (2) ≥ e− √ dR/λ. (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 R) 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.e.
E[KLCN,i,j ] ≤ EV (R)[KLCN,i,j ] (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 √ 3R, adding an error of O(e−2 √ 3R/λ). We then (2)
lower bound the distances of both points to both landmarks by the closest distance to a landmark r = min{‖xi − xl1‖2, ‖xi − xl2‖2, ‖xj − xl1‖2, ‖xj − xl2‖2} and use Lemma A to obtain
EV (R)[KLCN,i,j ] (1) = EV (R)[KLCN, 2 landmarks,i,j ] +O(e−2 √ 3R/λ)
(2) ≤ EV (R) [ 2e−2r/λ
1 + e−2R/λ
] +O(e−2 √ 3R/λ)
= 2EV (R)[e−2r/λ]
1 + e−2R/λ +O(e−2
√ 3R/λ).
(30)
Assuming Euclideanness in V (R) we obtain
EV (R)[e−2r/λ] = 1
V (R) ∫ R 0 e−2r/λ dV (r) dr dr = d R ∫ R 0 e−2r/λrd−1 dr
= d
(2R/λ)d (Γ(d)− Γ(d, 2R/λ))
(31)

heading: D PROOF OF THEOREM 2
text: Note that this theorem does not use probabilistic arguments but rather geometrically analyzes the maximum possible error. Ksp is correct for all pairs inside a cluster and 0 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 −Ksp‖∞ = e−D/λ (32) LCN is also correct for all pairs inside a cluster, so we again consider the closest possible pair xi, xj between clusters. We furthermore only consider the landmarks of the two concerned clusters, adding an error of O(e−D/λ). Hence, KLCN, 2 landmarks,ij = ( e−r/λ e−(r+D)/λ )( 1 e−(2r+D)/λ e−(2r+D)/λ 1 )−1( e−(r+D)/λ e−r/λ
) = 1
1− e−(4r+2D)/λ ( e−r/λ e−(r+D)/λ )( 1 −e−(2r+D)/λ −e−(2r+D)/λ 1 )( e−(r+D)/λ e−r/λ ) = 1
1− e−(4r+2D)/λ ( e−r/λ e−(r+D)/λ )(e−(r+D)/λ − e−(3r+D)/λ) e−r/λ − e−(3r+2D)/λ ) = 1
1− e−(4r+2D)/λ (e−(2r+D)/λ − e−(4r+D)/λ + e−(2r+D)/λ − e−(4r+3D)/λ)
= e−(2r+D)/λ
1− e−(4r+2D)/λ (2− e−2r/λ − e−(2r+2D)/λ)
= e−D/λe−2r/λ(2− e−2r/λ)−O(e−2D/λ) (33)
and thus max ‖K −KLCN‖∞ = e−D/λ(1− e−2r/λ(2− e−2r/λ) +O(e−D/λ)). (34)
For pure Nyström we need to consider the distances inside a cluster. In the worst case two points overlap, i.e. Kij = 1, and lie at the boundary of the cluster. Since r D we again only consider the landmarks in the concerned cluster, adding an error of O(e−D/λ). Because of symmetry we can optimize the worst-case distance from all landmarks by putting them on a (l− 1)-simplex centered on the cluster center. Since there are at most d landmarks in each cluster there is always one direction in which the worst-case points are r away from all landmarks. The circumradius of an (l − 1)-simplex with side length ∆ is √ l−1 2l ∆. Thus, the maximum distance to all landmarks is √ r2 + l−12l ∆
2. Using Lemma A we therefore obtain the Nyström approximation
KNys,ij = le−2 √ r2+ l−12l ∆ 2/λ
1 + (l − 1)e−∆/λ +O(e−D/λ) (35)

heading: E NOTE ON THEOREM 3
text: Lemmas C-F and and thus Theorem 1 by Altschuler et al. (2019) are also valid for Q outside the simplex so long as ‖Q‖1 = n and it 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 L22 in this theorem, which only changes the dependence on ρ.

heading: F NOTES ON THEOREM 4
text: To adapt Theorem 1 by Dvurechensky et al. (2018) to sparse matrices (i.e. matrices with some Kij = 0) we need to redefine
ν := min i,j {Kij |Kij > 0}, (36)
i.e. take the minimum only w.r.t. non-zero elements in their Lemma 1.

heading: G PROOF OF PROPOSITION 1
text: Theorem A (Danskin’s theorem). Consider a continuous function φ : Rk × Z → R, with the compact set Z ⊂ Rj . If φ(x, z) is convex in x for every z ∈ Z and φ(x, z) has a unique maximizer z̄, the derivative of
f(x) = max z∈Z φ(x, z) (37)
is given by the derivative at the maximizer, i.e.
∂f ∂x = ∂φ(x, z̄) ∂x . (38)
We start by deriving the derivatives of the distances. To show that the Sinkhorn distance fulfills the conditions for Danskin’s theorem we first identify x = C, z = P , and φ(C,P ) = −〈P ,C〉F + λH(P ). We next observe that the restrictions P1m = p and P T1n = q define a compact, convex set for P . Furthermore, φ is a continuous function and linear in C, i.e. both convex and concave for any finite P . Finally, φ(C,P ) is concave in P since 〈P ,C〉F is linear and λH(P ) is concave. Therefore the maximizer P̄ is unique and Danskin’s theorem applies to the Sinkhorn distance. Using
∂CNys,ij ∂Ukl = ∂ ∂Ukl
( −λ log(
∑ a UiaWaj)
) = −λδik
Wlj∑ aUiaWaj = −λδik Wlj KNys,ij , (39)
∂CNys,ij ∂Wkl = ∂ ∂Wkl
( −λ log(
∑ a UiaWaj)
) = −λδjl
Uik∑ aUiaWaj = −λδjl Uik KNys,ij , (40)
P̄Nys,ij KNys,ij = ∑ b P̄U,ibP̄W,bj∑ aUiaWaj = s̄it̄j ∑ bUibWbj∑ aUiaWaj = s̄it̄j ∑ bUibWbj∑ aUiaWaj = s̄it̄j (41)
and the chain rule we can calculate the derivative w.r.t. the cost matrix as
∂dλc ∂C = − ∂ ∂C
( −〈P̄ ,C〉F + λH(P̄ ) ) = P̄ , (42)
∂dλLCN,c ∂Ukl = ∑ i,j ∂CNys,ij ∂Ukl ∂dλLCN,c ∂CNys,ij = −λ ∑ i,j δikWlj P̄Nys,ij KNys,ij
= −λ ∑ i,j δikWlj s̄it̄j = −λs̄k ∑ j Wlj t̄j = ( −λs̄(Wt̄)T ) kl ,
(43)
∂dλLCN,c ∂Wkl = ∑ i,j ∂CNys,ij ∂Wkl ∂dλLCN,c ∂CNys,ij = −λ ∑ i,j δjlUik P̄Nys,ij KNys,ij
= −λ ∑ i,j δjlUiks̄it̄j = −λ (∑ i s̄iUik ) t̄l = ( −λ(s̄TU)T t̄T ) kl ,
(44)
and ∂dλLCN,c ∂ log Ksp and ∂dλLCN,c ∂ log KspNys directly follow from ∂d λ c ∂C . Theorem B (Implicit function theorem). Let f : Rn′ ×Rm′ → Rm′ be a continuously differentiable function with f(a, b) = 0. If its Jacobian matrix Jfi(x,y),yj = ∂fi ∂yj (a, b) is invertible, then there exists an open set a ∈ U ⊂ Rn′ on which there exists a unique continuously differentiable function g : U → Rm′ with g(a) = b and ∀x ∈ U : f(x, g(x)) = 0. Moreover,
∂gi ∂xj
(x) = − ( Jf(x,y),y(x, g(x)) )−1 i,:
( ∂
∂x f(x, g(x)) ) :,j
(45)
Next we derive the transport plan derivatives. To apply the implicit function theorem we identify x = C and y = P as flattened matrices. We will index these flat matrices via index pairs to simplify interpretation. We furthermore identify
f(C,P ) = ∂
∂P̄
( 〈P̄ ,C〉F − λH(P̄ ) ) = C + λ(log P̄ + 1). (46)
The minimizer P̄ cannot lie on the boundary of P ’s valid region since limp→0 ∂∂pp log p = limp→0 log p + 1 = −∞ and therefore P̄ij > 0. Hence, f(C, P̄ ) = 0 with P̄ (C) = arg minP 〈P ,C〉F − λH(P ) and we find that g(x) = P̄ (C). We furthermore obtain
Jfij(C,P ),Pab(C, P̄ (C)) = ∂
∂P̄ab
( Cij + λ(log P̄ij + 1) ) = δiaδjbλ/P̄ij , (47)
∂
∂Ckl fab(C, P̄ (C)) =
∂
∂Ckl
( Cab + λ(log P̄ab + 1) ) = δakδbl. (48)
Jfij(C,P ),Pab(C, P̄ (C)) is hence a diagonal matrix and invertible since λ/P̄ij > 0. We can thus use the implicit function theorem and obtain
∂P̄ij ∂Ckl = − ∑ a,b ( Jfij(C,P ),Pab(C, P̄ (C)) )−1 ∂ ∂Ckl fab(C, P̄ (C))
= − ∑ a,b δiaδjb 1 λ P̄ijδakδbl = − 1 λ P̄ijδikδjl.
(49)
To extend this result to LCN-OT we use
∂P̄U,ij ∂P̄Nys,ab = ∂ ∂P̄Nys,ab ∑ k P̄Nys,ikP̄ † W,kj = δiaP̄ † W,bj , (50)
∂P̄W,ij ∂P̄Nys,ab = ∂ ∂P̄Nys,ab ∑ k P̄ †U,ikP̄Nys,kj = δjbP̄ † U,ia (51)
and the chain rule to obtain
∂P̄U,ij ∂Ukl = ∑ a,b,c,d ∂CNys,cd ∂Ukl ∂P̄Nys,ab ∂CNys,cd ∂P̄U,ij ∂P̄Nys,ab
= ∑ a,b,c,d ( −λδck Wld KNys,cd )( − 1 λ P̄Nys,abδacδbd )( δiaP̄ † W,bj ) = ∑ b Wlb P̄Nys,ib KNys,ib P̄ †W,bjδik = δik ∑ b siWlbtbP̄ † W,bj = δiksi ∑ b P̄W,lbP̄ † W,bj
= δikδjlsi,
(52)
∂P̄W,ij ∂Ukl = ∑ a,b,c,d ∂CNys,cd ∂Ukl ∂P̄Nys,ab ∂CNys,cd ∂P̄W,ij ∂P̄Nys,ab
= ∑ a,b,c,d ( −λδck Wld KNys,cd )( − 1 λ P̄Nys,abδacδbd )( δjbP̄ † U,ia ) = Wlj
P̄Nys,kj KNys,kj P̄ †U,ik = WljtjP̄ † U,iksk = P̄ † U,ikskP̄W,lj ,
(53)
∂P̄U,ij ∂Wkl = ∑ a,b,c,d ∂CNys,cd ∂Wkl ∂P̄Nys,ab ∂CNys,cd ∂P̄U,ij ∂P̄Nys,ab
= ∑ a,b,c,d ( −λδdl Uck KNys,cd )( − 1 λ P̄Nys,abδacδbd )( δiaP̄ † W,bj ) = Uik
P̄Nys,il KNys,il P̄ †W,lj = siUiktlP̄ † W,lj = P̄U,iktlP̄ † W,lj ,
(54)
∂P̄W,ij ∂Wkl = ∑ a,b,c,d ∂CNys,cd ∂Wkl ∂P̄Nys,ab ∂CNys,cd ∂P̄W,ij ∂P̄Nys,ab
= ∑ a,b,c,d ( −λδdl Uck KNys,cd )( − 1 λ P̄Nys,abδacδbd )( δjbP̄ † U,ia ) = ∑ a Uak P̄Nys,aj KNys,aj P̄ †U,iaδjl = δjl ∑ a saUaktjP̄ † U,ia = δjltj ∑ a P̄ †U,iaP̄U,ak
= δikδjltj .
(55)
∂P̄ sp ∂ log Ksp and ∂P̄ spNys ∂ log KspNys directly follow from ∂P̄∂C .
We can calculate the pseudoinverses P̄ †U = (P̄ T U P̄U ) −1P̄U and P̄ † W = P̄ T W (P̄W P̄ T W ) −1 in time O((n + m)l2 + l3) since P̄U ∈ Rn×l and P̄W ∈ Rl×m. We do not fully instantiate the matrices required for backpropagation but instead save their decompositions, similar to the transport plan P̄Nys = P̄U P̄W . We can then compute backpropagation in time O((n + m)l2) by applying the sums over i and j in the right order. We thus obtain O((n + m)l2 + l3) overall runtime for backpropagation.

heading: H CHOOSING LSH NEIGHBORS AND NYSTRÖM LANDMARKS
text: We focus on two LSH methods for obtaining near neighbors. Cross-polytope LSH (Andoni et al., 2015) uses a random projection matrix R ∈ Rd×b/2 with the number of hash buckets b, and then decides on the hash bucket via h(x) = arg max([xTR ‖ −xTR]), where ‖ denotes concatenation. K-means LSH computes k-means and uses the clusters as hash buckets.
We further improve the sampling probabilities of cross-polytope LSH via the AND-OR construction. In this scheme we calculate B · r hash functions, divided into B sets (hash bands) of r hash functions
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, 2006).
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 LCNOT. For embedding alignment we use cross-polytope LSH 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 a better alignment between LSH samples and Nyström. We convert the cosine similarity to a distance via
dcos =
√ 1− x T p xq
‖xp‖2‖xq‖2 (Berg et al., 1984) to use k-means with dot product similarity. Note that
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 a GPU via PyTorch (Paszke et al., 2019) and PyTorch Scatter (Fey & Lenssen, 2019). To avoid over- and underflows we use log-stabilization throughout, i.e. we save all values in log-space and compute all matrix-vector products and additions via the log-sum-exp trick log ∑ i e xi = maxj xj + log( ∑ i e xi−maxj xj ). 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 or not calculating A−1 and instead solving the equation B = AX forX . We furthermore precomputeW = A−1V to avoid unnecessary operations.
We use 3 layers and an embedding size HN = 32 for GTN. The MLPs use a single hidden layer, biases and LeakyReLU non-linearities. The single-head MLP uses an output size of HN, match = HN and a hidden embedding size of 4HN, i.e. the same as the concatenated node embedding, and the multi-head MLP uses a hidden embedding size of HN. To stabilize initial training we scale the node embeddings by d̄
n̄ √ HN, match directly before calculating OT. d̄ denotes the average graph distance in the
training set, n̄ the average number of nodes per graph, and HN, match the matching embedding size, i.e. 32 for single-head and 128 for multi-head OT.

heading: J GRAPH DATASET GENERATION AND EXPERIMENTAL DETAILS
text: The dataset statistics are summarized in Table 5. Each dataset contains the distances between all graph pairs in each split, i.e. 10 296 and 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-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 4 and 1/2n̄ and 3/2n̄ total nodes, where n̄ 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̄/20 times by randomly adding, type editing, or removing nodes
and edges. Editing nodes and edges is 4x and adding/deleting edges 3x as likely as adding/deleting nodes. Most of these numbers were chosen arbitrarily, aiming to achieve a 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 we choose an edit cost of 1 for changing a node or edge type and 2 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 on 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 to 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 8 for the final values. We use the originally published result of SimGNN on Linux and thus don’t provide its hyperparameters. GTN uses 500 Sinkhorn iterations. We obtain the final entropy regularization parameter from λbase via λ = λbase d̄ n̄ 1 logn , where d̄ denotes the average graph distance and n̄ the average number of nodes per graph in the training set. The factor d̄/n̄ serves to estimate the embedding distance scale and 1/ log n counteracts the entropy scaling with n log n. Note that the entropy regularization parameter was small, but always far from 0, which shows that entropy regularization actually has a positive
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 is multiple times slower due to its iterative refinement scheme. Note that similarity matrix computation time (K) primarily depends on the LSH/Nyström method, not the OT approximation.
N = 10000 N = 20000 N = 50000 K OT K OT K OT
Full Sinkhorn 8 2950 29 11 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 LCN-Sinkhorn 101 116 242 205 642 624
0 100 200
Neighbors + landmarks
0
100
200
300
R un
tim e
(m s)
Multsc. OT Nys. Skh. Sparse Skh. LCN-Skh.
Figure 5: Runtime scales linearly with the number of neighbors/landmarks for all relevant Sinkhorn approximation methods.
100 1000 10000 Avg. graph size
10
100
1000
10000
Ti m
e pe
re po
ch (s
)
Full LCN
Figure 6: Log-log runtime per epoch for GTN with full Sinkhorn and LCN-Sinkhorn. LCNSinkhorn scales almost linearly with graph size while sustaining similar accuracy.
effect on learning. On the pref. att. 200 dataset we use no L2 regularization, λbase = 0.5, and a batch size of 200. For pref. att. 2k we use λbase = 2 and a batch size of 20 for full Sinkhorn and 100 for LCN-OT. For pref. att. 20k we use λbase = 50 and a batch size of 4. λbase scales with graph size due to normalization of the PM kernel.
For LCN-OT we use 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 is higher and strongly varies per iteration due to the unbalanced nature of hierarchical k-means. This increase in neighbors and landmarks and PyTorch’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.

heading: K RUNTIMES
text: Table 9 compares the runtime of the full Sinkhorn distance with different approximation methods using 40 neighbors/landmarks. We separate the computation of 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 is 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
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 as many multiplications for the same number of neighbors/landmarks. In either case, both Table 9 and Fig. 5 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 LCN-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 to sustain good accuracy. Full Sinkhorn went out of memory for the largest graphs.

heading: L DISTANCE APPROXIMATION
text: Figs. 7 and 8 show that for the chosen λ = 0.05 sparse Sinkhorn offers the best trade-off between computational budget and distance approximation, with 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λc can be negative due to the entropy offset. This picture changes as we increase the regularization. For higher regularizations LCN-Sinkhorn is the most precise at constant computational budget (number of neighbors/landmarks), as shown in Fig. 9. Note that the crossover points in 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 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.

Review summaries:
The 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.