{'fraction_non_alphanumeric': 0.0374899220639613, 'fraction_numerical': 0.014736182029920272, 'mean_word_length': 4.292012325195544, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Approximation of Wasserstein distance with Transshipment

An algorithm for approximating the p-Wasserstein distance between histograms defined on unstructured discrete grids is presented. It is based on the computation of a barycenter constrained to be supported on a low dimensional subspace, which corresponds to a transshipment problem. A multi-scale strategy is also considered. The method provides sparse transport matrices and can be applied to large scale and non structured data.


Introduction
The computation of optimal transportation between two discrete normalized mesaures µ x and µ y defined from x ∈ Ω ⊂ R d to [0; 1] remains a challenging problem when an accurate discretization of the domain Ω is considered. Optimal transportation requires to define a ground distance between points x, y ∈ R d of the domain Ω. This ground metric is then used to measure how much it costs to move µ x (x) to µ y (y). Ground distances ||x − y|| p are here considered for p ≥ 1, which leads to the p−Wasserstein distance between µ x and µ y (see [43,34,31] for more detailed introductions). Such distances give robust metrics in retrieval applications for 1 ≤ p < 2 [33,29,30,20]. The underlying sparse transport matrix is also a useful tool for interpolation and transfer purposes [32,16,40].
Computing Wasserstein distances The computation of Wasserstein distances is only explicit for d = 1. A standard approach to estimate Wasserstein distances when d > 1 consists in precomputing a cost matrix ||x − y|| p for x, y ∈ Ω and then estimating the whole transport matrix which dimension grows quadratically with n = |Ω|. Linear programming or transportation simplex [24,18] can be applied to estimate p−Wasserstein distance for p ≥ 1 but they are limited in practice to low dimensional problems, i.e. small values of n, for complexity and storage issues. The implementation provided in [10] nevertheless allows to tackle problems of interesting dimensions, i.e. Ω discretized with more than n = 10 4 points. By exploiting the sparsity of the transport map, multi-scale strategy [26] or grid refinement [36] can deal with larger problems by solving iteratively sparse low dimensional problems with linear programming. More efficient algorithms can be considered in the specific case of 1−Wasserstein distances [22].

Approximation of Wasserstein distance
For large scale problems, approximated sparse transport matrices and Wasserstein distances can be obtained by considering successive one dimensional problems with the so-called sliced Wasserstein distance [32,9]. The entropic regularization of the transport map proposed in [13] is another relevant way to deal with problems of high dimension. Given a regularization parameter γ > 0, it approximates the true Wasserstein distance (that corresponds in this setting to γ = 0) with the well-known Sinkhorn algorithm. The estimated transport matrices are nevertheless dense and should be truncated with care for interpolation purposes. When data are discretized on an uniform grid, the estimation of the distance can be obtained through iterative convolutions [40] with a kernel K p = exp(−||x − y|| p /γ). This leads to very fast algorithm for the p−Wasserstein distance as the Kernel is separable when considering the L p norm ||.|| p , and only convolutions and storage of d one-dimensional kernels are necessary. For small values of γ, numerical instabilities nevertheless arise and dedicated attention must be given to the implementation by considering for instance decaying values of γ, stabilization in log domain or multi-resolution approaches [37]. The numerical convergence is also reduced with low values of γ, and overrelaxation [35,41] or greedy coordinate descents [3] approaches have been proposed to tackle this issue. When the kernel K p is not separable, it must be carefully truncated to zero to store sparse matrices. If not considering more complex and adaptive truncations [37], it limits the possible amplitude of the transport, which may be annoying in case of large displacements. Other regularizations can then be of interest [15,8,39], namely when it is suitable to recover sparse transport matrices.
However, as mentioned before, this kind of techniques only leads to efficient implementations for densities discretized on structured grids. Moreover, if large displacements are involved in the data, the Wasserstein distance can not be accurately approximated, since numerical instabilities arises with ||x − y|| p /γ when γ goes to 0.

Wasserstein barycenters and Transshipment
The p−Wasserstein barycenter, as introduced in [1], can be used to approximate the p−Wasserstein distance between measures involving a large number n of dirac masses. As proposed in [44] for clustering problems, the discrete barycenter between 2 discrete measures can be parameterized with a weighted sum of κ << n diracs to obtain a low dimensional problem that corresponds to the transshipment problem of n resources with κ intermediate locations. The sum of the distances between each data and the barycenter then gives an approximation of the effective distance, as it has been proposed for 1−Wasserstein distances [5]. Statistical properties of such a method have been later studied in [17,28], where it has been underlined the robustness of this low rank regularization of the transport matrix to data outliers. These ideas have been extended to discrete approximation of barycenter between continuous measures in [11].
Content The use of low dimensional barycenters is the point of view adopted in this note to propose fast approximation of Wasserstein distances involving sparse transport matrices. As in [44], at each iteration of the presented algorithm, the κ locations of the barycenter are updated and a linear program of dimension 2κn is solved. Compared to the dimension n 2 for the classical Wasserstein distance, the overall complexity of the algorithm becomes attractive for high dimensional data. In [44] and [17,28], this low dimensional barycenter problem is respectively solved with the Alternating Direction Method of Multipliers or Sinkhorn iterations [14]. The barycenter problem between two discrete densities is here seen as a transshipment problem and solved up to numerical accuracy with an efficient network simplex graph algorithm based on the work of [10]. Then κ transportation subproblems are solved in parallel to recover a sparse transport matrix. Following [37,26,23], this leads to the design of a multi-scale barycenter scheme, to iteratively refine the transport matrix and the associated approximation of the Wasserstein distance.
The theoretical computer science community has recently provided improved bounds for graph based algorithms solving approximate transportation problems up to additive or 1 + multiplicative errors. This namely includes the optimal transportation problem [7], its entropic regularization [2] or the transshipment problem [6]. As for the sliced method [32], the presented approach does not have such theoretical guarantees but very good performances are observed in practice. More precisely, approximate p−Wasserstein distances between cloud points of 10 5 elements are obtained in a few minutes without involving prohibitive memory storage issues. The proposed empirical algorithm can also be directly applied to non structured data. Thanks to the multi-scale refinement approach, a sparse transport map is provided which can be of interest for interpolation purposes.
Outline Section 2 details how to approximate the Wasserstein distance when computing the barycenter between two discrete measures. The strategy of [44] is recalled in the general context of p−Wasserstein distance. A multi-scale algorithm for recovering sparse transport matrices is finally presented. The performances of the algorithm are discussed in section 3 through extensive experiments realized on the benchmark [38]. Numerical results show that this whole empirical process is efficient, namely when one of the two input data is spatially regular.

Approximate Wasserstein distance from barycenter estimation
Let µ x and µ y be two discrete measures defined on They have normalized positive weights vectors w x ∈ S m and w y ∈ S n , where S n is the simplex of size n defined as be the ground cost matrix over Ω defined as c xy ij = x i − y j p , which corresponds to the power p of the distance related to a given norm on R d . Then the p−Wasserstein distance between discrete measures µ x and µ y is with the set of admissible transport matrices and where 1 n the vector full of ones in R n . This problem can be efficiently solved with linear programming. It can also been formulated through a directed graph containing m + n nodes and mn vertices. The final transport matrix γ is very sparse in practice (at most m + n + 1 non null entries) but the involved complexity and memory storage scale with the product of data dimensions nm. For latter purpose, Algorithm 1 details the function estimating the distance W p p (µ x , µ y ).

Interpolation and barycenters
Let γ xy be an optimal transport matrix solution of (1) and µ t be the following interpolation between measures µ x and µ y for t ∈ [0; 1] : This interpolation is the discrete analogue [31] to the geodesic between µ x and µ y defined by the McCann's interpolation [25]. Discrete versions of some results in [34] can now be expressed.
Proposition 1 For p ≥ 1, γ xy a solution of (1), µ t defined in (3) and t ∈ [0; 1], the following relations hold: Proof. Observing that involves that the previous relations are in fact equalities.
Following [1], it can be shown that the mid interpolation µ 1/2 is solution of the p-Wasserstein barycenter problem between µ x and µ y with weights (1/2, 1/2). Proposition 2 For p ≥ 1 and the interpolation µ 1/2 defined in (3), it holds that and µ 1/2 is a solution of the p−Wasserstein barycenter problem: Proof. From (4), it can first be noticed that: Since W p is a distance and the function |.| p , is convex for p ≥ 1, it can next be observed that ∀µ: Combing relations (8) and (9) implies that µ 1/2 is a solution of the barycenter problem (7).
Existence (and uniqueness for p > 1) of Wasserstein barycenters have been deeply studied in [1,19]. As stated in the following proposition, the p−Wasserstein distance can be obtained for any p ≥ 1 through the resolution of the barycenter problem (7).

Barycenter computation
A discrete barycenter µ 1/2 between µ x and µ y can be obtained by solving (7) with the distance (1). As done in [14] with an additional entropic regularization, the mid barycenter µ 1/2 between µ x and µ y can be constrained to be supported on a set of κ dirac masses, i.e. µ 1/2 = κ k=1 w z k δ z k , with positions z k ∈ R d and weights w z ∈ S κ . The barycenter problem can be rewritten as: with the admissible set of matrices and the cost matrices c xz ik = x i − z k p and c yz jk = y j − z k p . Notice that the weight vector of the barycenter w z ∈ R κ is implicitly included in the set of constraints: (γ x ) 1 m = (γ y ) 1 n = w z . The positions z k here act as intermediate locations where the mass has to transit from x to y. Following [44] and as illustrated in Figure 1, the idea behind this modeling is to consider a limited number of transshipment locations κ to speed up the computation.
Optimal Transportation problem Optimal Transshipment problem The problem (12) is separately convex with respect its variables γ x , γ y and z k . For p > 1, the coupling terms are differentiable and alternate minimization over transport matrices (γ x , γ y ) and dirac positions z k converges [42] to a saddle point. Following (10), the Wasserstein distance W p (µ x , µ y ) can be approximated withW p (µ x , µ y ) obtained from transport matrices γ x and γ y solutions of the problem (12) as The alternate optimization steps for solving (12) are now detailed and the process is summed up in Algorithm 2.
Remark 1 The norm . 2 is usually taken as reference for the p−Wasserstein distance. From the equivalence of norms in finite dimensions and as can be done with entropic regularization to make the problem more tractable numerically, the L p norm x i −y j p p = d s=1 (x s i −y s j ) p is here considered for computing the p−Wasserstein distance. Get γ x and γ y by solving transshipment problem (15) under the constraints (16) 7: Setz k = z k 8: Update positions z k from γ x and γ y by following Sec. 2.2.1 9: until z k −z k / z k <

Update of positions
To update positions z k ∈ R d , the problem (12) is solved for fixed transport matrices γ x and γ y . For p ≥ 1, this leads to d convex problems, that can be solved in parallel for each space dimension s=1. . . d. Dimension indexes s are thus omitted in the following and the problem writes Different strategies are considered according to p. Some locations z k may become useless if i γ x ik = j γ y jk = 0. In this case the corresponding z k are removed and κ is decreased. When considering large scale problems and few transshipment locations κ, this almost never happens.
Case W 2 . For p = 2, there exists an explicit update formula of the dirac positions to find the unique minimizer of (13) with respect to z: This step acts like the cluster position update in a κ-mean algorithm. It realizes for each z k a weighted mean of the positions x i 's and y j 's according to γ x ik and γ y jk .
Case W p , p > 2. For p > 2 the problem (13) admits a unique minimizer and is twice differentiable. Newton's method can be considered to approximate the solution: where c xz ik = z k − x i and c yz jk = z k − y j , while ⊗ denotes the elementwise product between matrices and the power p − 2 is also element-wise.
Case W 1 . When p = 1, a global optimum of (13) can be obtained by taking each z k as a weighted median of the positions x i and y j with respect to the weights γ x ik and (γ x ) jk . This operation mainly requires to sort the value of the x i 's and y j 's along each dimension. Notice that alternate minimization on problem (12) may not converge when p = 1.
Case W p , p ∈]1; 2[. The problem does not admits a second derivative. An iterative scheme is then considered by decomposing |z k − x i | p = |z k − x i | p−1 |z k − x i | and solving successive weighted median problems between the positions x i and y j with the weights γ x ik |z k − x i | p−1 and γ y jk |z k − y j | p−1 .

Update of transport matrices through Transshipment
For fixed z k , problem (12) can be solved with classic linear programming optimization tools. The interesting point is that this barycenter problem is a transshipment problem with κ intermediate locations. It can therefore be formulated in terms of a directed graph with m + κ + n vertices (i.e. x i , z k and y j ) and (mκ + κn) edges e ik and e kj that correspond to the transport matrices γ x ik and γ y jk . The following cost function is then minimized under the set of constraintP(w x , w y ) that translates into: This problem can be efficiently solved with the network simplex algorithm [27]. An extension of the original non sparse implementation proposed in [10] is here considered.

Discussion
It is well known [1,4] that if µ x and µ y are respectively supported by m and n dirac masses, then their exists a barycenter supported by up to m + n + 1 dirac masses. Hence, the approximation of the barycenter from a set of κ < min(m, n) dirac masses seems to be a bad choice at first sight. Notice however that, as pointed out in [17], a barycenter supported by a low dimensional space allows the approximate distance to be more robust to data outliers. Next, contrary to entropic regularization of OT, the obtained approximate transport map is here sparse, which allows efficient storage and can be directly used for interpolation purposes without any complex post-processing like sharpening.  Finally, it is worth noting that the approximation depends on the regularity of the data. As underlined in [1], when computing the barycenter of a set of densities µ i in the continuous case, if one of the input data µ i is absolutely continuous with respect to the Lebesgue measure, so does the barycenter. Such observation gives an interesting insight of the experiments presented in this paper: when at least one of the two data µ x or µ y is smooth, then the obtained approximate Wasserstein distanceW p can be very close to the true one W p for small values of κ. This point is illustrated in Figure 2 with the comparison of relative errors between true and approximate distances obtained with increasing values of κ for different scenarios involving "random" or "smooth" 2D data of dimensions m = n = 10 4 . When at least one "regular" data is involved, then a small relative error (W p − W p )/W p < 10 −3 (i.e. 0.1%) is observed for κ/n = 0.015.

Transport refinement and Multi-scale approach
As illustrated in Figure 2, the barycenter approach is not sufficient to produce an accurate approximation of the p−Wasserstein distance between any data. First notice that the approximationW p of the Wasserstein distance can be easily improved using the following result that directly considers the transport matrix between µ x and µ y .
Proposition 4 Let µ z = κ k=1 w z k δ z k and γ x (resp. γ y ) be an optimal transport matrix from µ x (resp. µ y ) to µ z . Let alsoγ be the transport matrix between µ x and µ y defined asγ = γ x D(γ y ) , with the rescaling D −1 = diag(w z ) given by the weights w z k = i γ x ik = j γ y jk . Then the following relation holds  Proof. The matrixγ ij = k γ x ik γ y jk /w z k is an admissible transport matrix between µ x and µ y , since jγ ij = k γ x ik = w x i and iγ ij = k γ y jk = w y j . Following the proof of Minkowski's inequality: , then gives (17). Notice that the proof of the triangle inequality [12] is here also obtained, since From this proposition, the approximationŴ p (µ x , µ y ) corresponding to a barycenter µ 1/2 = κ k=1 w z k δ z k solution of (12) with transport matrices γ x and γ y is defined as For matching or interpolation purposes,γ = γ x (diag(w z )) −1 (γ y ) gives a sparse approximation of the optimal transport matrix. However, as illustrated in Figure 3, such approach maps all locations x i and y j that transit by z k , resulting in a poor block transport matrix for small values of κ. µ x µ y True transport κ = 16 κ = 32 κ = 128 Figure 3: Illustration of the block transport matricesγ estimated from data µ x and µ y through transshipment for increasing values of κ.
In order to get a sparser and more accurate approximation of the optimal transport matrix and associated distance, it is necessary to "untie" the links between locations passing through z k . To do so, a solution is to consider κ optimal transportation sub-problems, by refining the transport of the mass transshipped through z k . The barycenter approach then acts as a clustering and for each intermediate location z k , the p−Wasserstein distance between the following partial discrete densities is computed: where the number of active dimensions are expected to be reduced: m k = #{γ x ik > 0} < m and n k = #{γ y jk > 0} < n. The multi-scale approach is then performed as follows. The Wasserstein distance W p (µ k x , µ k y ) is estimated exactly with network simplex [10] if the sum of m k and n k is small enough. Otherwise the barycenter approach is recursively applied to the subproblem. From numerical experiments, the threshold n k +m k < 2000 has been chosen to reach the best compromise between numerical accuracy and computational cost. The whole process is illustrated in Figure 4 and detailed in Algorithm 3.

3:
InitializeŴ = 0,γ = 0 m×n 4: for k = 1 to κ do Can be done in parallel 5: Set m k = #I k , n k = #J k 7: if m k + n k < T then 9:  A more detailed presentation is then given in Figure 5, with the mean relative errors over intra and inter classes experiments in the case of 128 × 128 images. As can be observed, significant errors are obtained when unstructured random data (classes 4 and 10) are involved. In all other cases, the relative errors are very low.  Computational cost In Figure 6, the running time for computing an approximate 2−Wasserstein distance, with Algorithm 3 and κ = 16, is compared with the C ++ network simplex implementation and its sparse multi-threaded extension proposed in [10]. The running times become asymptotically very interesting with the multi-threaded extension, but due to memory storage, these methods can not handle dimensions n larger than 3.10 4 on the considered computer. As transsshipment involves problems of size nκ, it can be applied to data containing more dirac masses and thus deal with one additional order of magnitude (n = 2.10 5 ). With the proposed full C ++ implementation of the transshipment problem, the κ sub-problems are solved successively. The provided Matlab interface calling the C ++ code through a simple parfor loop with 4 workers is thus much faster. Optimal transshipment matrices being dense for small values of κ, it is counter-productive to consider sparse optimized implementation in the multi-scale framework.

Conclusion
This paper presents an empirical method for approximating Wasserstein distances. It is based on existing concepts also used in parallel works [17,5,28]. The contribution is to provide an efficient multi-scale implementation able to deal with unstructured point clouds while providing sparse transport matrices. Numerical experiments demonstrate the accuracy of the computed approximate distance, while the involved computational cost are improved with respect to the literature. As a perspective, it would be of interest to add constraints encouraging an homogeneous repartition of the number of points transiting by each location z k of the barycenter, or at least a more uniform distribution of the barycenter weights w z k . The κ multi-scale sub-problems that can be solved in parallel would thus have similar dimensions, and theoretical guarantees on the overall running time could be given. Parallel Approximate W p Figure 6: Comparison of running times for computing 2−Wassersetin distances for different values of n: C ++ network simplex [10], its multi-threaded extension and the proposed multi-scale approximation (matlab C ++ mex), that is the only one being able to handle high scale problems.