Efficient estimates of optimal transport via low-dimensional embeddings Conference Submissions
Year: 2020

Abstract:
Optimal transport distances (OT) have been widely used in recent work in Machine Learning as ways to compare probability distributions. These are costly to compute when the data lives in high dimension. Recent work aims specifically at reducing this cost by computing OT using low-rank projections of the data (seen as discrete measures) (Paty & Cuturi, 2019). We extend this approach and show that one can approximate OT distances by using more general families of maps provided they are 1-Lipschitz. The best estimate is obtained by maximising OT over the given family. As OT calculations are done after mapping data to a lower dimensional space, our method scales well with the original data dimension. We demonstrate the idea with neural networks. We use Sinkhorn Divergences (SD) to approximate OT distances as they are differentiable and allow for gradient-based optimisation. We illustrate on synthetic data how our technique preserves accuracy and displays a low sensitivity of computational costs to the data dimension.

heading: 1 Introduction
text: Optimal Transport metrics (Kantorovich, 1960) or Wasserstein distances, have emerged successfully in the field of machine learning, as outlined in the review by Peyré et al. (2017). They provide machinery to lift distances on X to distances over probability distributions in P(X ). They have found multiple applications in machine learning: domain adaptation Courty et al. (2017), density estimation (Bassetti et al., 2006) and generative networks (Genevay et al., 2017; Patrini et al., 2018). However, it is prohibitively expensive to compute OT between distributions with support in a high-dimensional space and might not even be practically possible as the sample complexity can grow exponentially as shown by Dudley (1969). Similarly, work by Weed et al. (2019) showed a theoretical improvement when the support of distributions is found in a low-dimensional space. Furthermore, picking the ground metric that one should use is not obvious when using high-dimensional data. One of the earlier ideas from Santambrogio (2015) showed that OT projections in a 1-D space may be sufficient enough to extract geometric information from high dimensional data. This further prompted Kolouri et al. (2018) to use this method to build generative models, namely the Sliced Wasserstein Autoencoder. Following a similar approach Paty & Cuturi (2019) and Muzellec & Cuturi (2019) project the measures into a linear subspace E of low-dimension k that maximizes the transport cost and show how this can be used in applications of color transfer and domain adaptation. This can be seen as an extension to earlier work by Cuturi & Doucet (2014) whereby the cost function is parameterized.
One of the fundamental innovations that made OT appealing to the machine learning com-
munity was the seminal paper by Cuturi (2013) that introduced the idea of entropic regularization of OT distances and the Sinkhorn algorithm. Since then, regularized OT has been successfully used as a loss function to construct generative models such as GANs (Genevay et al., 2017) or RBMs (Montavon et al., 2015) and computing Barycenters (Cuturi & Doucet, 2014; Claici et al., 2018). More recently, the new class of Sinkhorn Divergences was shown by Feydy et al. (2018) to have good geometric properties, and interpolate between Maximum Mean Discrepancies (MMD) and OT.
Building on this previous work, we introduce a general framework for approximating highdimensional OT using low-dimensional projections f by finding the subspace with the worst OT cost, i.e. the one maximizing the ground cost on the low-dimensional space. By taking a general family of parameterizable fφs that are 1-Lipschitz, we show that our method generates a pseudo-metric and is computationally efficient and robust. We start the paper in §2 with background on optimal transport and pseudo-metrics. In §3 we define the theoretical framework for approximating OT distances and show how both linear (Paty & Cuturi, 2019) and non-linear projections can be seen as a special instance of our framework. In §4 we present an efficient algorithm for computing OT distances using Sinkhorn Divergences and fφs that are 1-Lipschitz under the L2 norm. We conclude in §5 with experiments illustrating the efficiency and robustness of our method.

heading: 2 Preliminaries
text: We start with a brief reminder of the basic notions needed for the rest of the paper. Let X be a set equipped with a map dX : X × X → R≥0 with non-negative real values. The pair (X , dX ) is said to be a metric space and dX is said to be a metric on X if it satisfies the usual properties:
• dX (x, y) = 0 if and only if x = y • dX (x, y) = dX (y, x) • dX (x, z) ≤ dX (x, y) + dX (y, z)
If dX verifies the above except for the only if condition, it is called a pseudo-metric, and (X , dX ) is said to be a pseudo-metric space. For a pseudo-metric, it may be that dX (x, y) = 0 while x 6= y. We write dX ≤ d′X if for all x, y dX (x, y) ≤ d′X (x, y). It is easy to see that: 1) “≤” is a partial order on pseudo-metrics over X, 2) “≤” induces a complete lattice structure on the set of pseudo-metrics over X , where 3) suprema are computed pointwise (but not infima). Consider X , Y, two metric spaces equipped with respective metrics dX , dY . A map f from X to Y is said to be α-Lipschitz continuous if dY(f(x), f(x′)) ≤ αdX (x, x′). A 1-Lipschitz map is also called non-expansive. Given a map f from X to Y one defines the pullback of dY along f as:
f̂(dY)(x, x′) = dY(f(x), f(x′)) (1)
It is easily seen that: 1) f̂(dY) is a pseudo-metric on X , 2) f̂(dY) is a metric iff f is injective, 3) f̂(dY) ≤ dX iff f is non-expansive, 4) f̂(dY) is the least pseudo-metric on the set X such that f is non-expansive from (X , f̂(dY)) to (X , dY). Thereafter, we assume that all metric spaces considered are complete and separable, i.e. have a dense countable subset.
Let (X , dX ) be a (complete separable) metric space. Let ΣX be the σ-algebra generated by the open sets of X (aka the Borelian subsets). We write P(X ) for the set of probability distributions on (X ,ΣX). Given a measurable map f : X → Y, and µ ∈P(X) one defines the push-forward of µ along f as:
f#(µ)(B) = µ(f−1(B)) (2) for B ∈ ΣY . It is easily seen that f#(µ) is a probability measure on (Y,ΣY ) Given µ in P(X ), ν in P(Y), a coupling of µ and ν is a probability measure γ over X × Y such that for all A in ΣX , B in ΣY , γ(A×X ) = µ(A), and γ(X ×B) = ν(B). Equivalently, µ = π0#(π), and ν = π1#(π) for π0, π1 the respective projections. We write Γ(µ, ν) for the set of couplings of µ and ν. There are several ways to lift a given metric structure on dX to one on P(X ). We will be specifically interested in metrics on P(X ) derived from optimal transport problems. The p-Wasserstein metric with p ∈ [1,∞) is defined by:
Wp(dX )(µ, ν)p = inf γ∈Γ(µ,ν) ∫ X×X dpX dγ (3)
Villani (2008) establishes that if dX is (pseudo-) metric so isWp(dX ). The natural ‘Dirac’ embedding of X into P(X ) is isometric (there is only one coupling). The idea behind the definition is that dpX is used as a measure of the cost of transporting units of mass in X , while a coupling γ specifies how to transport the µ distribution to the ν one. One can therefore compute the mean transportation cost under γ, and pick the optimal γ. Hence the name optimal transport. In most of the paper, we are concerned with the case X = Rd+ for some large d with a metric structure dX given by the Euclidean norm, and we wish to compute the W2 metric between distributions with finite support. Since OT metrics are costly to compute in high dimension, to estimate these efficiently, and mitigate the impact of dimension, we will use a well-chosen family of fs to push the data along a map with a low dimensional co-domain Y also equipped with the Euclidean metric. The reduction maps may be linear or not. They have to be non-expansive to guarantee that the associated pull-back metrics are always below the Euclidean one, and therefore we provide a lower estimate of W2(d2).

heading: 3 Approximate OT with General Projections - GPW
text: With the ingredients from the above section in place, we can now construct a general framework for approximating Wasserstein-like metrics by low-dimensional mappings of X . We write simply W instead of Wp as the value of p plays no role in the development. Pick two metric spaces (X , dX ), (Y, dY), and a family S = (fφ : X → Y; φ ∈ S) of mappings from X to Y. Define a map from P(X )×P(X ) to non-negative reals as follows:
dS(µ, ν) = sup S W (dY)(fφ#(µ), fφ#(ν)) (4)
Equivalently and more concisely dS can be defined as: dS(µ, ν) = sup
φ W (f̂φ(dY))(µ, ν) (5)
It is easily seen that:
1. the two definitions are equivalent 2. dS is a pseudo-metric on P(X ) 3. dS is a metric (not just a pseudo one) if the family fφ jointly separates points in X , and 4. if the fφs are non-expansive from (X , dX ) to (Y, dY), then dS ≤W (dX )
The second point follows readily from the second definition. Each f̂φ(dY) is a pseudo-metric on X obtained by pulling back dY (see preceding section), hence, so is W (f̂φ(dY)) on P(X ), and therefore dS being the supremum of this family (in the lattice of pseudo-metrics over X ) is itself a pseudo-metric. The first definition is important because it allows one to perform the OT computation in the target space where it will be cheaper. Thus we have derived from S a pseudo-metric dS on the space of probability measures P(X ). We assume from now on that mappings in S are non-expansive. By point 4. above, we know that dS is bounded above by W (dX ). We call dS the generalized projected Wasserstein metric (GPW) associated to S. In good cases, it is both cheaper to compute and a good estimate.

heading: 3.1 SRW as an instance of GPW
text: In Paty & Cuturi (2019), the authors propose to estimate W2 metrics by projecting the ambient Euclidean X into k-dimensional linear Euclidean subspaces. Specifically, their derived metric on P(X), written Sk, can be defined as (Paty & Cuturi, 2019, Th. 1, Eq. 4):
S2k(µ, ν) = sup Ω W 22 (dY)(Ω 1/2 # (µ),Ω 1/2 # (ν)) (6)
where: 1) dY is the Euclidean metric on Y, 2) Ω contains all positive semi-definite matrices of trace k (and therefore admitting a well-defined square root) with associated semi-metric smaller than dX . We recognise a particular case of our framework where the family of mappings is given by the linear mappings √ Ω : Rd = X → Y = Rk under the constraints above. In particular, all mappings used are linear. The authors can complement the general properties of the approach with a specific explicit bound on the error and show that S2k ≤ W 22 (dX ) ≤ (d/k)S2k. In the general case, there is no upper bound available, and one has only the lower one.

heading: 3.2 Non-linear embeddings for approximating Wasserstein distances
text: Using the same Euclidean metric spaces, X = Rd, Y = Rk, we observe that our framework does not restrict us to use linear functions as mappings. One could use a family of mapping given by a neural network (fφ : X → Y;φ ∈ S) where φ ranges over network weights. However, not any φ is correct. Indeed, by point 4) in the list of properties of dS , we need fφs to be non-expansive. Ideally, we could pick S to be the set of all weights such that fφ is non-expansive. There are two problems one needs to solve in order to reduce the idea to actual tractable computations. First, one needs an efficient gradient-based search to look for the weights φ which maximise supSW (dY)(fφ#(µ), fφ#(ν)) (see 4). Second, as the gradient update may take the current fφ out of the non-expansive maps, one needs to project back efficiently in the space of non-expansive. Both problems already have solutions which are going to re-use. For the first point, we will use Sinkhorn Divergence (SD) (Genevay et al., 2017). Recent work Feydy et al. (2018) shows that SD, which one can think as a regularised version of W , is a sound choice as a loss function in machine
learning. It can approximate W closely and without bias (Genevay et al., 2017), has better sample complexity (Genevay et al., 2019), as well as quadratic computation time. Most importantly, it is fully differentiable. For the second problem, one can ‘Lipshify’ the linear layers of the network by dividing their (operator) norm after each update. We will use linear layers with Euclidean metrics, and this will need to estimate the spectral radius of each layer. The same could be done with linear layers using a mixture of L1, L2 and L∞ metrics. In fact computing the L1 → L1 operator norm for linear layers is an exact operation, as opposed to using the spectral norm for L2 → L2 case where we approximate using the power method. Note that the power method can only approximate the L2 norm and gradient ascent methods used in the maximization phase are stochastic making our approximation susceptible to more variables. However, it is extremely efficient since it requires computation of optimal transport distances only in the low-dimensional space. We can see this as a trade-off between exactness and efficiency.

heading: 4 Computational details
text: In this section we propose Algorithm 1 for stochastically estimating dS between two measures with finite support where the class of mappings S is as defined above. Note that this algorithm can further be used during the training of a discriminator as part of a generative network with an optimal transport objective, similar to Genevay et al. (2017). The Sinkhorn Divergence alternative for dS now uses Sinkhorn divergences as a proxy for OT (compare with equation 4):
SDφ, (µ, ν) = W (dY)(fφ#(µ), fφ#(ν))− 1 2W (dY)(fφ#(µ), fφ#(µ))− 1 2W (dY)(fφ#(ν), fφ#(ν))
(7)
where W is the well-known Sinkhorn regularized OT problem (Cuturi, 2013). The nonparameterized version of the divergence has been shown by Feydy et al. (2018) to be an unbiased estimator of W (µ, ν) and converges to the true OT distance when = 0. Their paper also constructs an effective numerical scheme for computing the gradients of the Sinkhorn divergence on GPU, without having to back-propagate through the Sinkhorn iterations, by using autodifferentiation and the detach methods available in PyTorch (Paszke et al., 2019). Moreover, work by Schmitzer (2019) devised an -scaling scheme to trade-off between guaranteed convergence and speed. This gives us further control over how fast the algorithm is. It is important to note that the minimization computation happens in the low-dimensional space, differently from the approach in Paty & Cuturi (2019), which makes our algorithm scale better with dimension, as seen in §5.
Feydy et al. (2018) established that the gradient of 7 w.r.t to the input measures µ, ν is given by the dual optimal potentials. Since we are pushing the measures through a differentiable function fφ, we can do the maximization step via a stochastic gradient ascent method such as SGD or ADAM (Kingma & Ba, 2014). Finally, after each iteration, we project back into the space of 1-Lipschitz functions fφ. For domain-codomain L2 ←→ L2 the Lipschitz constant of a fully connected layer is given by the spectral norm of the weights, which can be approximated in a few iterations of the power method. Since non-linear activation functions such as ReLU are 1-Lipschitz, in order to project back into the space of constraints we suggest to normalize each layer’s weights with the spectral norm, i.e. for layer i we have φi := φi/||φi||. Previous work done by Neyshabur et al. (2017) as well as Yoshida & Miyato (2017) and Miyato et al. (2018) showed that with smaller magnitude weights, the model can better generalize and improve the quality of generated samples when used on a discriminator in a GAN. We note that if we let fφ to be a 1-Layer fully connected network with no activation, the optimization we perform is very similar with the optimization done by Paty & Cuturi (2019). The space of 1-Lipschitz functions we are optimizing over is larger and
our method is stochastic, but we are able to recover very similar results at convergence. Moreover, our method applies to situations where the data lives in a non-linear manifold that an fφ such as a neural network is able to model. Comparing different numerical properties of the Subspace Robust Wasserstein distances in 6 with our Generalized Projected Wasserstein Distances is the focus of the next section.
Algorithm 1 Ground metric parameterization through φ Input: Measures µ = ∑n i δxiai and ν = ∑n j δyj bj , fφ : Rd → Rk 2-Layer network with dimen-
sions (d, 20, k) and 1-Lipschitz, optimizer ADAM , power method iterations λ, SDφ, unbiased Sinkhorn Divergence. Output: fφ, SDφ, Initialize: lr, ,λ, fφ ∼ N (0, 10), Objective← SD (blur = 2, p = 2, debias = True) for t→ 1, . . . ,maxiter do L← −SDφ, (fφ#µ, fφ#ν) (pushforward through fφ and evaluate SD in lower space) gradφ ← Autodiff ( L) (maximization step with autodiff)
φ← φ+ADAM(gradφ) (gradient step with SGD and scheduler) φ← Projλ1−Lip(φ) (projection into 1-Lipschitz space of functions)
end for

heading: 5 Experiments
text: We consider similar experiments as presented in Forrow et al. (2019) and Paty & Cuturi (2019) and show the mean estimation of SD2φ,k(µ, ν) for different values of k, as well as robustness to noise. We also show how close the distance generated by the linear projector from Paty & Cuturi (2019) is to our distance and highlight the trade-off in terms of computation time with increasing number of dimensions. In order to illustrate our method, we construct two empirical distributions µ̂, ν̂ by taking samples from two independent measures µ = N (0,Σ1) and ν = N (0,Σ2) that live in a 10 dimensional space. Similarly to Paty & Cuturi (2019) we construct the covariance matrices Σ1,Σ2 such that they are of rank 5, i.e. the support of the distributions is given by a 5 dimensional linear subspace. Throughout our experiments we fix fφ to be a 2-layer neural network with a hidden layer of 16 units, activation function ReLU and output of dimension k. We initialize the weights from N (0, 10) and use a standard ADAM optimizer with a decaying cyclic learning rate (Smith, 2017) bounded by [0.1, 1.0]. Decreasing and increasing the learning rate via a scheduler allows us to not fall into local optima. The batch size for the algorithm is set to n = 500, which is the same number of samples that make up the two measures. Besides the neural network variables, we set the regularization strength small enough, to = 0.001, and the scaling to -scaling = 0.95 such that we can accurately estimate the true optimal transport distance, but not spend too much computational time during the Sinkhorn iterates.

heading: 5.1 10-D Gaussian Data OT estimation using SDφ,k
text: This leaves us with three variables of interest during the computation of SDφ,k, namely k, d, λ (latent dimension, input dimension, power method iterations). The power method iterations plays an important role during the projection step, as for a small number of iterations, there is a chance
2 3 4 5 6 7 8 9 10 k
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
M ea
n No
rm al
ize d
di st
an ce
s
Noisy SD2( , ) W2( , ) SD2( , ) W2( , ) Noisy S2k ( , ) W2( , ) S2k ( , ) W2( , )
Figure 2: Mean normalized distances with and without noise for SD2φ(µ, ν) and S2k(µ, ν) as a function of latent dimension k. The shaded area shows the standard deviation over 20 runs.
of breaking the constraint. At the same time, running the algorithm for too long is computationally expensive. In Figure 1 we used λ = 5 power iterations and show the values of SD2k,φ after running 1 for 500 iterations. We compare them to the true OT distance for various levels of k and observe that even with a small number of power iterations, the estimation approaches the true value as k increases. Furthermore, we see that for k = 5 and k = 7 the algorithm converges after 200 steps.
Using 20 power iterations, we show how the approximation behaves in the presence of noise as a function of the latent space k. We add Gaussian noise in the form of N (0, I) to µ̂, ν̂ and show in Figure 2 the comparison between no noise and noise for both SRW distances defined in 6 and GPW in 4. We observe that SD2φ,k behaves similarly to S2k in the presence of noise.

heading: 5.2 Computation time
text: In Figure. 8 of Paty & Cuturi (2019) they note that their method when using Sinkhorn iterates is quadratic in dimension because of the eigen-decomposition of the displacement matrix. Fundamentally different, we are always optimizing in the embedded space, making the computation of the Sinkhorn iterates linear with dimension. Note that there is the extra computation involved with pushing the measures through the neural network and backpropagating as well as the projection step that depends on the power iteration method. In order to run this experiment we set λ = 5 and generate µ̂, ν̂ by changing dimension d but leaving the rank of Σ1,Σ2 equal to 5. The latent space is fixed to k = 5. In Figure 3 we plot the normalized distances using the two approaches as a function of dimension and see that the gap gets bigger with increasing dimensions, but it is stable. In Figure 4 we plot the log of the relative computation time, taking the d = 10 as a benchmark in both cases. We see that the time to compute SD2φ is linear in dimension and is significantly lower than its counterpart S2k as we increase the number of dimensions. This can be traced back to Algorithm. 1 and Algorithm. 2 of Paty & Cuturi (2019) where at each iteration step, the computation of OT distances in the data space is prohibitively expensive.
0 250 500 750 1000 1250 1500 1750 2000 Dimension
100
101
Re la
tiv e
co m
pu ta
tio n
tim e
(lo g
sc al e) SD2( , ) S2k ( , )
Figure 4: Mean relative computation time (log scale) comparison between the two distances. The shaded area shows shows the standard deviation over 20 runs.

heading: 6 Conclusion
text: In this paper we presented a new framework for approximating optimal transport distances using a wide family of embedding functions that are 1-Lipschitz. We showed how linear projectors can be considered as a special case of such functions and proceeded to define neural networks as another class of embeddings. We showed how we can use existing tools to build an efficient algorithm that is robust and constant in the dimension of the data. Future work includes showing the approximation is valid for datasets where the support of distributions lies in a low-dimensional non-linear manifold, where we hypothesize that linear projects would fail. Other work includes experimenting with different operator norms such as L1 or Linf for the linear layers and the approximation of W1. An extension of the projection step in 1 to convolutional layers would allow us to experiment with real datasets such as CIFAR-10 and learn a discriminator in an adversarial way with SDk,φ as a loss function. This can be used to show that the data naturally clusters in the embedding space.

Review summaries:
This paper uses the fact that the Wasserstein distance is decreasing under 1-Lipschitz mappings from the ambient space to a feature space in order to propose more robust (to dimensionality ?) estimation of the Wasserstein distance between two probability distributions. A neural network with some sort of weight renormalization is used to produce 1-Lipschitz embeddings. The authors then maximise over the parametrized maps.