Year: 2021

Abstract:
Optimal Transport (OT) metrics allow for defining discrepancies between two probability measures. Wasserstein distance is for longer the celebrated OT-distance frequently-used in the literature, which seeks probability distributions to be supported on the same metric space. Because of its high computational complexity, several approximate Wasserstein distances have been proposed based on entropy regularization or on slicing, and one-dimensional Wassserstein computation. In this paper, we propose a novel extension of Wasserstein distance to compare two incomparable distributions, that hinges on the idea of distributional slicing, embeddings, and on computing the closed-form Wassertein distance between the sliced distributions. We provide a theoretical analysis of this new divergence, called heterogeneous Wasserstein discrepancy (HWD), and we show that it preserves several interesting properties including rotation-invariance. We show that the embeddings involved in HWD can be efficiently learned. Finally, we provide a large set of experiments illustrating the behavior of HWD as a divergence in the context of generative modeling and in query framework.

heading: 1 INTRODUCTION
text: Optimal Transport-based data analysis has recently found widespread interest in machine learning community, since its significant usefulness to achieve many tasks arising from designing loss functions in supervised learning (Frogner et al., 2015), unsupervised learning (Arjovsky et al., 2017), text classification (Kusner et al., 2015), domain adaptation (Courty et al., 2017), generative models (Arjovsky et al., 2017; Salimans et al., 2018), computer vision (Bonneel et al., 2011; Solomon et al., 2015) among many more applications (Kolouri et al., 2017; Peyré & Cuturi, 2019). Optimal Transport (OT) attempts to match real-world entities through computing distances between distributions, and for that it exploits prior geometric knowledge on the base spaces in which the distributions are valued. Computing OT distance equals to finding the most cost-efficiency way to transport mass from source distribution to target distribution, and it is often referred to as the Monge-Kantorovich or Wasserstein distance (Monge, 1781; Kantorovich, 1942; Villani, 2009).
Matching distributions using Wasserstein distance relies on the assumption that their base spaces must be the same, or that at least a meaningful pairwise distance between the supports of these distributions can be computed. A variant of Wasserstein distance dealing with heterogeneous distributions and overcoming the lack of intrinsic correspondence between their base spaces is Gromov-Wasserstein (GW) distance (Sturm, 2006; Mémoli, 2011). GW distance allows to learn an optimal transport-like plan by measuring how the similarity distances between pairs of supports within each ground space are closed. It is increasingly finding applications for learning problems in shape matching (Mémoli, 2011), graph partitioning and matching (Xu et al., 2019), matching of vocabulary sets between different languages (Alvarez-Melis & Jaakkola, 2018), generative models (Bunne et al., 2019), or matching weighted networks (Chowdhury & Mémoli, 2018). Due to the heterogeneity of the distributions, GW distance uses only the relational aspects in each domain, such as the pairwise relationships to compare the two distributions. As a consequence, the main disadvantage of GW distance is its computational cost as the associated optimization problem is a non-convex quadratic program (Peyré & Cuturi, 2019), and as few as thousand samples can be computationally challenging. Based on the approach of regularized OT (Cuturi, 2013), in which an entropic penalty is added to the original objective function defining the Wasserstein OT problem, Peyré et al. (2016) propose an entropic version called entropic GW discrepancy, that leads to approximate GW distance. Another
approach for scaling up the GW distance is Sliced Gromov-Wasserstein (SGW) discrepancy (Vayer et al., 2019), which leverages on random projections on 1D and on a closed-form solution of the 1D-Gromov-Wasserstein.
In this paper, we take a different approach for measuring the discrepancy between two heterogeneous distributions. Unlike GW distance that compares pairwise distances of elements from each distribution, we consider a method that embeds the metric measure spaces into a one-dimensional space and computes a Wasserstein distance between the two 1D-projected distributions. The key element of our approach is to learn two mappings that transform vectors from the unit-sphere of a latent space to the unit-sphere of the metric space underlying the two distributions of interest, see Figure 1. In a nutshell, we learn to transform a random direction, sampled under an optimal (learned) distribution (optimality being made clear later), from a d-dimensional space to a random direction into the desired spaces. This approach has the benefit of avoiding an ad-hoc padding strategy (completion of 0 of the smaller dimension distributions to fit the high-dimensional one) as in SGW method (Vayer et al., 2019). Another relevant feature of our approach is that the two resulting 1D distributions are now compared through Wasserstein distance. This point, in conjunction, with other key aspect of the method, will lead to a relevant discrepancy between two distributions, called heterogeneous Wasserstein discrepancy (HWD). Although we lose some properties of a distance, we show that HWD is rotation-invariant, that it is robust enough to be considered as a loss for learning generative models between heterogeneous spaces. We also establish that HWD boils down to the recent distributional sliced Wasserstein distance (Nguyen et al., 2020) if the two distributions live in the same space and if some mild constraints are imposed on the mappings.
In summary, our contributions are as follows:
• we propose HWD, a novel slicing-based discrepancy for comparing two distributions living in different spaces. Our chosen formulation is based on comparing 1D random-projected versions of the two distributions using a Wasserstein distance;
• The projection operations are materialized by optimally mapping from one common space to the two spaces of interest. We provide a theoretical analysis of the resulting discrepancy and exhibit its relevant properties;
• Since the discrepancy involves several mappings that need to be optimized, we depict an alternate optimization algorithm for learning them;
• Numerically, we validate the benefits of HWD in terms of comparison between heterogeneous distributions. We show that it can be used as a loss for generative models or shape objects retrieval with better performance and robustness than SGW on those tasks.

heading: 2 BACKGROUND OF OT DISTANCES
text: For the reader’s convenience, we provide here a brief review of the notations and definitions, that will be frequently used throughout the paper. We start by introducing Wasserstein and Gromov-Wasserstein distances with their sliced versions SW and SGW, where we consider these distances in the specific case of Euclidean base spaces (Rp, ‖ · ‖) and (Rq, ‖ · ‖). We denote P(X ) and P(Y) the respective sets of probability measures whose supports are contained on compact sets X ⊆ Rp and Y ⊆ Rq. For r ≥ 1, we denote Pr(X ) the subset of measures in P(X ) with finite r-th moment (r ≥ 1), i.e., Pr(X ) = { η ∈ P(X ) :´
X ‖x‖ rdη(x) < ∞
} . For µ ∈ P(X ) and ν ∈ P(Y), we write Π(µ, ν) ⊂ P(X × Y)
for the collection of joint probability distributions with marginals µ and ν, known as couplings, Π(µ, ν) = { γ ∈P(X × Y) : ∀A ⊂ X , B ⊂ Y, γ(A× Y) = µ(A), γ(X ×B) = ν(B) } .

heading: 2.1 OT DISTANCES FOR HOMOGENEOUS DOMAINS
text: We here assume that the distributions µ and ν lie in the same base space, for instance p = q. Taking this into account, we can define the Wasserstein distance and its sliced variant.
Wasserstein distance The r-th Wasserstein distance is defined on Pr(X ) by
Wr(µ, ν) = (
inf γ∈Π(µ,ν) ˆ X×Y ‖x− y‖rdγ(x, y) ) 1 r . (1)
The quantityWr(µ, ν) describes the least amount effort to transform one distribution µ into another one ν. Since the cost distance used between sample supports is the Euclidean one, the infimum in (1) is attained (Villani, 2009), and any probability γ which realizes the minimum is called an optimal transport plan. In a finite discrete setting, Problem (1) can be formulated as a linear program, that is challenging to solve algorithmically as its computational cost is of order O(n5/2 log n) (Lee & Sidford, 2014), where n is the number of sample supports.
Contrastingly, for the 1D case (i.e. p = 1) of continuous probability measures, the r-th Wasserstein distance has a closed-form solution (Rachev & Rüschendorf, 1998), namely, Wr(µ, ν) = ( ´ 1
0 |F−1µ (u)− F−1ν (u)|rdt) 1 r where F−1µ and F −1 ν are the quantile functions of µ and ν. For empirical distributions, the 1D-Wasserstein distance is simply calculated by sorting the supports of the distributions on the real line, resulting to a complexity of order O(n log n). This nice computational property motivates the use of sliced-Wasserstein (SW) distance (Rabin et al., 2012; Bonneel et al., 2015), where one calculates an (infinity) of 1D-Wasserstein distances between linear projection pushforwards of distributions in question and then computes their average.
To precisely define SW distance, we consider the following notation. Let Sp−1 := {u ∈ Rp : ‖u‖ = 1} be the unit sphere in p dimension in `2-norm, and for any vector θ in Sp−1, we define Pθ the orthogonal projection onto the real line Rθ = {αθ : α ∈ R}, that is Pθ(x) = 〈θ, x〉, where 〈·, ·〉 stands for the Euclidean inner-product. Let µθ = Pθ#µ the measure on the real line called pushforward of µ by Pθ, that is µθ(A) = µ(P−1θ (A)) for all Borel set A ⊆ R. We may now define the SW distance.
Sliced Wasserstein distance The r-th order sliced Wasserstein distance between two probability distributions µ, ν ∈Pr(X ) is given by
SWr(µ, ν) = ( 1 Ap ˆ Sp−1 Wrr (µθ, νθ)dθ ) 1 r , (2)
where Ap is the area of the surface of Sp−1, i.e., Ap = 2π p/2
Γ(p/2) with Γ : R → R, the Gamma function given as Γ(u) = ´∞ 0 tu−1e−tdt. Thanks to its computational benefits and its valid metric property (Bonnotte, 2013), the SW distance has recently been used for OT-based deep generative modeling (Kolouri et al., 2019; Deshpande et al., 2019; Wu et al., 2019). Note that the normalized integral in (2) can be seen as the expectation for θ ∼ σp−1, the uniform surface measure on Sp−1, that is SWr(µ, ν) = (Eθ∼σp−1 [Wrr (µθ, νθ)]) 1 r . Therefore, the SW distance can be easily approximated via a Monte Carlo sampling scheme by drawing uniform random samples from Sp−1:
SWrr(µ, ν) ≈ 1K ∑K k=1Wrr (µθk , νθk) where θ1, . . . , θK
i.i.d.∼ σp−1 and K is the number of random projections.

heading: 2.2 OT DISTANCES FOR HETEROGENEOUS DOMAINS
text: To get benefit from the advantages of OT in many machine learning applications involving heterogeneous and incomparable domains (p 6= q), the Gromov-Wasserstein distance (Mémoli, 2011) stands for the basic OT distance dealing with this setting.
Gromov-Wasserstein distance The r-th Gromov-Wasserstein distance between two probability distributions µ ∈Pr(X ) and ν ∈Pr(Y) is defined by
GWr(µ, ν) = inf γ∈Π(µ,ν)
Jr(γ) def. =
1
2 ( ¨ X 2×Y2 |‖x− x′‖ − ‖y − y′‖|rdγ(x, y)dγ(x′, y′) ) 1 r . (3)
Note that GWr(µ, ν) is a valid metric endowing the collection of all isomorphism classes metric measure spaces of Pr(X )×Pr(Y), see Theorem 5 in (Mémoli, 2011). The GW distance learns an optimal transport-like plan which transports samples from a source metric space X into a target metric space Y , by measuring how the similarity distances between pairs of samples within each space are close. Furthermore, GW distance enjoys several geometric properties, particularly translation and rotation invariance. However, its major bottleneck consists in an expensive computational cost, since problem (3) is non-convex and quadratic. A remedy to such a heavy computational burden lies in an entropic regularized GW discrepancy (Peyré et al., 2016), using Sinkhorn iterations algorithm (Cuturi, 2013). This latter needs a large regularization parameter to guarantee a fast computation, which, unfortunately, entails a poor approximation of the true GW distance value. Another approach to scale up the computation of GW distance is sliced-GW discrepancy (Vayer et al., 2019). The definition of SGW shows 1D-GW distances between projected pushforward of an artifact zero padding of µ or ν distribution. We detail this representation in the following paragraph.
Sliced Gromov-Wasserstein discrepancy Assume that p < q and let ∆ be an artifact zero padding from X onto Y , i.e. ∆(x) = (x1, . . . , xp, 0, . . . , 0) ∈ Rq. The r-th order sliced Gromov-Wasserstein discrepancy between two probability distributions µ ∈Pr(X ) and ν ∈Pr(Y) is given by
SGW∆,r(µ, ν) = ( Eθ∼σq−1 [ GWrr((∆#µ)θ, νθ) ]) 1r . (4)
It is worthy to note that SGW∆,r is depending on the ad-hoc operator ∆, hence the rotation invariance is lost. Vayer et al. (2019) propose a variant of SGW that does not depend on the choice of ∆, called Rotation Invariant SGW (RI-SGW) for p = q, defined as the minimizer of SGW∆,r over the Stiefel manifold, see (Vayer et al., 2019, Equation 6). In this work, we are interested in calculating an OT-based discrepancy between distributions over distinct domains using the slicing technique. Our approach is different from the SGW one in many points, specifically (and most importantly) we use a 1D-Wasserstein distance between the projected pushforward distributions and not a 1D-GW distance. In the next section, we detail the setup of our approach.

heading: 3 HETEROGENEOUS WASSERSTEIN DISCREPANCY
text: Despite the computational benefit of sliced-OT variant discrepancies, they have an unavoidable bottleneck corresponding to an intractable computation of the expectation with respect to uniform distribution of projections. Furthermore, the Monte Carlo sampling scheme can often generate an overwhelming number of irrelevant directions; hence, the larger number of sample projections, the more accurate approximation of sliced-OT values. Recently, Nguyen et al. (2020) have proposed the distributional-SW distance allowing to find an optimal distribution over an expansion area of informative directions. This performs the projection efficiently by choosing an optimal number of important random projections needed to capture the structure of distributions. Our approach for comparing distributions in heterogeneous domains follows a distributional slicing technique combined with OT metric measure embedding (Alaya et al., 2020).
Let us first introduce additional notations. Fix d ≥ 1 and consider two nonlinear mappings φ : Sd−1 → Sp−1 and ψ : Sd−1 → Sq−1. For any constants Cφ, Cψ > 0, we define the following probability measure sets: MCφ = { Θ ∈ P(Sd−1) : Eθ,θ′∼Θ[|〈φ(θ), φ(θ′)〉|] ≤ Cφ } and
MCψ = { Θ ∈ P(Sd−1) : Eθ,θ′∼Θ[|〈ψ(θ), ψ(θ′)〉|] ≤ Cψ } . We say that Cφ, Cψ are (φ, ψ)admissible constants if the intersection sets MCφ ∩ MCψ is not empty. We hereafter denote µφ,θ = Pφ(θ)#µ and νψ,θ = Pψ(θ)#ν the pushforwards of µ and ν by projections over unit sphere Pφ(θ) and Pψ(θ), respectively.
Informal presentation While the distributions µ and ν are valued in different spaces, X ⊂ Rp and Y ⊂ Rq, any projected distributions will live in real line, enabling the computation of 1D-Wasserstein distance (Figure 1, right). In order to generate random 1D projections in each of the spaces, we map a common random projection distribution from Sd−1 into each of the projection spaces Sp−1 and Sq−1, through the mappings φ and ψ (see Figure 1, left). Hence, the main components of the heterogeneous Wasserstein discrepancy will be the distribution Θ ∈ P(Sd−1), and the two embeddings φ and ψ which will be wisely chosen. The resulting directions φ(θ) and ψ(θ) form the projections Pφ(θ) and Pψ(θ) (see Figure 1, center) used to compute several 1D-Wasserstein distances.

heading: 3.1 DEFINITION AND PROPERTIES
text: Herein we state the formulation of the proposed discrepancy and exhibit its main theoretical properties.
Definition 1 The heterogeneous Wasserstein discrepancy (HWD) of order r ≥ 1 between µ ∈ Pr(X ) and ν ∈Pr(Y) reads as
HWDr(µ, ν) = inf φ,ψ sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr (µφ,θ, νψ,θ) ]) 1r . (5)
HWD belongs to a family of projected OT works (Paty & Cuturi, 2019; Rowland et al., 2019; Lin et al., 2021) with a particularity for seeking nonlinear projections minimizing a sliced-OT variant. HWD further inherits the distributional slicing benefit by finding an optimal probability measure Θ of slices on the unit sphere Sd−1 coupled with an optimum couple (φ, ψ) of embeddings. Note that this optimal Θ verifies the double conditions Eθ,θ′∼Θ[| cos(φ(θ), φ(θ′))|] ≤ Cφ and Eθ,θ′∼Θ[| cos(ψ(θ), ψ(θ′))|] ≤ Cψ. This gives that Cφ, Cψ ≤ 1, hence the sets MCφ and MCψ belong to M1 = {Θ ∈P(Sd−1)} the set of all probability measures of the unit sphere Sd−1. It is worthy to note that for small regularizing (φ, ψ)-admissible constants, the measure Θ is forced to distribute more weights to directions that are far from each other in terms of their angles (Nguyen et al., 2020).
Now, in order to guarantee the existence of (φ, ψ)-admissible constants, we assume that the couple(φ, ψ)-embeddings are approximately angle preserving.
Assumption 1 (Approximately angle preserving property) For any couple (φ, ψ)-embeddings , assume that there exists two non-negative constants Lφ and Lψ, such that the following holds
|〈φ(θ), φ(θ′)〉| ≤ Lφ|〈θ, θ′〉| and |〈ψ(θ), ψ(θ′)〉| ≤ Lψ|〈θ, θ′〉|, for all θ, θ′ ∈ Sd−1.
In Proposition 1, we deliver lower bounds of the regularizing (φ, ψ)-admissible constants Cφ and Cψ, depending on the dimension of the latent space d and on the levels (Lφ, Lψ) of approximately angle preserving property. These bounds ensure the non-emptiness of the sets MCφ and MCψ .
Proposition 1 Let Assumption 1 hold and consider regularizing (φ, ψ)-admissible constants such that Cφ ≥ LφΓ(d/2)√πΓ((d+1)/2) and Cψ ≥ LψΓ(d/2)√ πΓ((d+1)/2) . Then the sets MCφ and MCψ contain the uniform
measure σd−1 and σ̄ = ∑d k=1 1 dδθk , where {θ1, . . . , θd} forms any orthonormal basis in R
d. Note that by Gautschi’s inequality (Gautschi, 1959) for the Gamma function, we have that Cφ ≥ Lφd and Cψ ≥ Lψd .
Proof of Proposition 1 is presented in Appendix A.2. Together the admissible constants, the levels of angle preserving property, and the dimension d of the latent space form the hyperparameters set of
HWD problem. For settings of large d, the admissible constants could take smaller values, that force the measure Θ to focus on far-angle directions. However, for smaller d, we may lose the control on the distributional part, the set MC tends to M1 the entire set of probability measure on Sd−1, hence it boils down on a standard slicing approach that needs an expensive number of projections to get an accurate approximation. Next, we give a set of interesting theoretical properties characterizing HWD.
Proposition 2 For any r ≥ 1, HWD satisfies the following properties:
(i) HWDr(µ, ν) is finite, that isHWDr(µ, ν) ≤ 2 r−1 r (Mr(µ) +Mr(ν)) where Mr(·) is the r-th moment of the given distribution, i.e. Mr(µ) = ( ´ X ‖x‖ rdµ(x) ) 1 r .
(ii) HWDr(µ, ν) is non-negative, symmetric and verifiesHWDr(µ, µ) = 0.
(iii) HWDr(µ, ν) has a discrepancy equivalence given by(1 d ) 1 r inf φ,ψ max θ∈Sd−1 Wr(µφ,θ, νψ,θ) ≤ HWDr(µ, ν) ≤ inf φ,ψ max θ∈Sd−1 Wr(µφ,θ, νψ,θ).
(iv) For p = q, HWD is upper bounded by the distributional sliced Wasserstein distance.
(v) HWDr is rotation invariant, namely, HWDr(R#µ,Q#ν) = HWDr(µ, ν), for any R ∈ Op = {R ∈ Rp×p : R>R = Ip} and Q ∈ Oq = {Q ∈ Rq×q : Q>Q = Iq}, the orthogonal group of rotations of order p and q, respectively.
(vi) Let Tα and Tβ be the translations from Rp into Rp and from Rq into Rq with vectors α and β, respectively. ThenHWDr(Tα#µ, Tβ#ν) ≤ 2r−1 ( HWDr(µ, ν) + ‖α‖+ ‖β‖) ) .
Proof of Proposition 2 is given in Apprendix A.1. From property (i) HWD is finite provided that the distributions in question have a finite r-th moments. Note that the supremum over the probability measure sets MCφ and MCψ guarantees the propertyHWDr(µ, µ) = 0. For p = q, if the infimum over the couple (φ, ψ)-embedding in (iii) is realized in the identity mappings, then HWD verifies a metric equivalence with respect to the max-sliced Wasserstein distance (Deshpande et al., 2019). The property (v) highlights a rotation invariance of HWD, which is well verified by the GW distance.

heading: 3.2 ALGORITHM
text: Computing HWD requires a resolution of an optimization problem as given in (5). In what follows, we propose an algorithm for computing an approximation of this discrepancy based on samples X = {xi}ni=1 from µ and samples Y = {yj}mj=1 from ν. At first, let us note that we have min-max optimization to solve; the minimization occuring over the embeddings φ and ψ and the maximization over the distributions on the unit-sphere. This maximization problem is challenging due to both the constraints and because we optimize over distributions. Similarly to Nguyen et al. (2020), we approximate the problem by replacing the constraints with regularization terms and by replacing the optimization over distributions by an optimization over a push-forward of the uniform probability measure σd−1 by a Borel measurable function f : Sd−1 → Sd−1. Hence, assuming that we have drawn from a uniform distribution K directions {θk}Kk=1, the numerical approximation of HWD is obtained by solving the following problem:
min φ,ψ max f
{ L1 def. = ( 1 K K∑ k=1 Wrr ( X>φ[f(θk)], Y >ψ[f(θk)] ))1/r}
+min f λC { L2 def. = ∑ k,k′ φ[f(θk)] >φ[f(θk′)] + ∑ k,k′ ψ[f(θk)] >ψ[f(θk′)] } +min f λa { L3 def. = ∑ k,,k′ ( φ[f(θk)] >φ[f(θk′)]− θ>k θk′ )2 + ∑ k,k′ ( ψ[f(θk)] >ψ[f(θk′)]− θ>k θk′ )2}
where the first term in the optimization is related to the sliced Wasserstein, the second term is related to the regularization term associated to Eθ,θ′∼Θ[|〈φ(θ), φ(θ′)〉|], and Eθ,θ′∼Θ[|〈ψ(θ), ψ(θ′)〉|], and the third term is the angle-preserving regularization term. Note that the min term with respect to f is due to the fact that we want those regularizers to be small. Two hyperparameters λC and λa
Algorithm 1 COMPUTING HETEROGENEOUS WASSERSTEIN DISCREPANCY (SEE (5)) 1: Input: Source and target samples: (X,µ) and (Y, ν); order r; the set of random direction {θk}Kk=1; T number of global iterations; N number of iterations for each of the alternate scheme;
2: Output: HWD 3: function COMPUTE LOSSES(X , Y , φ[f(θk)] , ψ[f(θk)]) 4: compute the average of 1D-Wasserstein L1 between X>φ[f(θk)] and Y >ψ[f(θk)] 5: compute the inner-product penalty L2 6: compute the angle-preserving penalty L3 7: end function 8: for t = 1, · · · , T do 9: fix φ and ψ
10: for i = 1, · · · , N do 11: compute φ[f(θk)] and ψ[f(θk)] 12: L1, L2, L3 ← Compute Losses (X , Y , φ[f(θk)] , ψ[f(θk)])) 13: L = −L1 + L2 + L3 14: f ← f − γk∇L 15: end for 16: fix f 17: for i = 1, · · · , N do 18: compute φ[f(θk)] and ψ[f(θk)] 19: L1, L2, L3 ← Compute Losses (X,Y , φ[f(θk)] , ψ[f(θk)]) 20: L = L1 + L2 + L3 21: φ← φ− γk∇L 22: ψ ← ψ − γk∇L 23: end for 24: end for 25: HWD ← compute the average over {θk}Kk=1 of closed-form 1D Wasserstein between
X>φ[f(θk)] and Y >ψ[f(θk)] 26: Return: HWD
control the impact of these two regularization terms. In practice, φ, ψ and f are parametrized as deep neural networks and the min-max problem is solved by an alternating optimization scheme : (a) optimizing over f with ψ and φ fixed then (b) optimizing over ψ and φ with f fixed. Some details of the algorithms are provided in Algorithm 1.
Regarding computational complexity, if we assume that the mapping f, φ, ψ are already trained, that we have K projections, and that φ[f(θk)], ψ[f(θk)] are precomputed, then the computation of HWD (line 25 of Algorithm 1) is in O(K(n log n + np + nq)), where n is the number of samples in X and Y . When taking into account the full optimization process, then the complexity depends on the number of times we compute the full objective function we are optimizing. Each evaluation requires the computation of the sum in L1 which is O(K(n log n + np + nq)) and the two regularization terms L2 and L3 require both O(K2(p+ q + 2d)). Note that in terms of computational complexity, SGW is O(Kn log n) whereas HWD is O(TNKn log n), with T ×N being the global number of objective function evaluations. Hence, complexity is in favor of SGW. However, one should note that in practice, because we optimize over the distribution of the random projections, we usually need less slices than SGW and thus depending on the problem, TNK can be of the same magnitude than the number of slices involved in SGW (similar findings have been highlighted for Sliced Wasserstein distance (Nguyen et al., 2020)).

heading: 4 NUMERICAL EXPERIMENTS
text: In this section, we analyze HWD, exhibit its rotation-invariant property, and compare its performance with SGW in a generative model context.
Translation and Rotation We have used two simple datasets for showing the behavior of HWD with respect to translation and rotation. For translation, we consider two 2D Gaussian distributions one being fixed, the other with varying mean. For rotation, we use two 2D spirals from Scikit-Learn
Figure 2: Examples of distance computation between (left) two-translating Gaussian distributions. (right) two spirals.
4 3
2 1
0 1
2 3
4 4 3
2 1
0 1
2 3
4
x3
4
3
2
1
0
1
2
3
4
target
3 2 1 0 1 2 3 x1
3
2
1
0
1
2
3
x2
target
Figure 3: Examples of target distributions for our generative models (left) 3D 4-mode. (right) 2D 5-mode.
library, one being fixed and the other being rotated from 0 to π/2. For these two cases, we have drawn 500 samples, used 100 random directions for SGW and RI-SGW. For our HWD, we have used only 10 slices and T = 50, N = 5 iterations for each of the alternate optimization. The results we obtain are depicted in Figure 2. From the first panel, we remark that both SGW and RI-SGW are indeed insensitive to translation while HWD captures this translation, which is also verified by property (vi) in Proposition 2. For the spiral problem, as expected HWD and RI-SGW are indeed rotation-invariant while SGW is not.
Generative models For checking whether our distribution discrepancy behaves appropriately, we have used it as a loss function in a generative model. Our task here is to build a model able to generate a distribution defined on a space having a different dimensionality from the target distribution space. As such, we have considered the same toy problems as in Bunne et al. (2019) and investigated two situations: generating 2D distributions from 3D data and the other way around. The 3D target distribution is a Gaussian mixture model with four modes while the 2D ones are 5-mode. Our generative model is composed of a fully-connected neural network with ReLU activation functions. We have considered 3000 samples in the target distributions and batch size of 300. For both HWD and SGW, we have run the algorithm for 30000 iterations with an Adam optimizer, stepsize of 0.001 and default β parameters. For the 4-mode problem, the generator is a MLP with 2 layers while for the 5-mode, as the problem is more complex, it has 3 layers. In each case, we have 256 units on the first layer and then 128. For the hyperparameters, we have set λC = 1 and λa = 5 or λa = 50 depending on the problem. Note that for SGW, we have also added a `2-norm regularizer on the output of the generator in order to avoid them to drift (see Figures 8 and 9 in Appendix C), as the loss is translation-invariant. Examples of generated distributions are depicted in Figure 4. We remark that our HWD is able to produce visually correct distributions whereas SGW struggles in generating the 4 modes and its 3D 5-mode is squeezed on its third dimension.
Scalability We consider the non-rigid shape world dataset (Bronstein et al., 2006) which consists of 148 three-dimensional shapes from 12 classes. We draw randomly n ∈ {100, 250, 500, 1000, 1500, 2000} vertices {xi ∈ R3}ni=1 on each shape and use them to measure the similarity between a pair of shapes {xi ∈ R3}ni=1 and {yj ∈ R3}nj=1. Figure 5 reports the average time to compute on a single core such a similarity for 100 pairs of shapes using respectively GW, SGW and HWD. As expected GW exhibits a slow behavior while the computational burden of HWD is on par with SGW.
Classification under various transformations This experiment, whose details are provided in Appendix B.2, aims to evaluate the robustness of SGW and HWD (the computationally efficient methods) to different transformations in terms of classification accuracy. To that purpose we employ the Shape Retrieval Contest (SHREC’2010) correspondence dataset, see Bronstein et al. (2010). It includes high resolution (10K-50K) triangular meshes. The shapes are of 3 classes (see Figure 10 in Appendix C) with 9 different transformations and the null shape (no transformation). Each transformation is applied up to five strength levels (weak to strong). Along with the null shape, we consider all strengths of the "isometry", "topology", "scale", "shotnoise" transformations leading to 63 samples. We perform a 1- NN classification. Obtained performances over 10 runs are depicted
in Figure 6. They highlight the ability of HWD to be robust to perturbations. HWD achieves slightly better mean classification accuracy than SGW with a competitive computation time (see Figure 5). Notice that GW and RISGW are unable to run under reasonable time-budget constraint.

heading: 5 CONCLUSION
text: We introduce in this paper HWD a novel OT-based discrepancy between distributions lying in different spaces. It takes computational benefits from distributional slicing technique, which amounts to find an optimal number of random projections needed to capture the structure of data distributions. Another feature of this discrepancy consists in projecting the distributions in question through a learning of embeddings enjoying the same latent space. We showed a nice geometrical property verified by the proposed discrepancy, specifically a rotation-invariance. We illustrated through extensive experiments the applicability of this discrepancy on generative modeling and shape objects retrieval. We argue that the implementation part faces the standard deep learning bottleneck of tuning the model’s hyperparameters. A future extension line of this work is to deliver theoretical guarantees regarding the regularizing parameters, both of distributional and angle preserving properties.

heading: A PROOFS
text: A.1 PROOF OF PROPOSITION 1
We use the following result:
Lemma 1 [Theorem 3 in (Nguyen et al., 2020)] For uniform measure σd−1 on the unit sphere Sd−1, we have ˆ
Sd−1×Sd−1 |〈θ, θ′〉|dσd−1(θ)dσd−1(θ′) = Γ(d/2)√ πΓ((d+ 1)/2) .
Hence,ˆ Sd−1×Sd−1 |〈φ(θ), φ(θ′)〉|dσd−1(θ)dσd−1(θ′) (Assumption 1) ≤ Lφ ˆ Sd−1×Sd−1 |〈θ, θ′〉|dσd−1(θ)dσd−1(θ′)
(Lemma 1) ≤ LφΓ(d/2)√
πΓ((d+ 1)/2) .
Therefore, as long as the (φ, ψ)-admissible constants Cφ ≥ LφΓ(d/2)√πΓ((d+1)/2) and Cψ ≥ LψΓ(d/2)√ πΓ((d+1)/2) , we have σd−1 ∈MCφ ∩MCψ . Now using a Gautschi’s inequality (Gautschi, 1959) for the Gamma function, it yields that Γ(d/2)√
πΓ((d+1)/2) ≥ 1√ π(d+1)/2 ≥ 1/d. Let σ̄ = ∑d l=1 1 dδθl , where {θ1, . . . , θd}
form an orthonormal basis in Rd. We then have Eθ,θ′∼σ̄ [ |〈φ(θ), φ(θ′)〉| ] = ∑ 1≤k,l≤d (1 d )2|〈φ(θk), φ(θ′l)〉| (Assumption 1)≤ Lφ ∑ 1≤k,l≤d (1 d )2|〈θk, θ′l〉| = Lφd . Therefore we get the lower bounds for the (φ, ψ)-admissible constants Cφ and Cψ given in Proposition 1, that guarantee σd−1, σ̄ ∈MCφ ∩MCψ .
A.2 PROOF OF PROPOSITION 2
Let us first state the two following lemmas: Lemma 2 writes an integration result using push-forward measures; it relates integrals with respect to a measure η and its push-forward under a measurable map f : X → Y. Lemma 3 proves that the admissible set of couplings between the embedded measures are exactly the embedded of the admissible couplings between the original measures.
Lemma 2 [See Lerner (2014) p. 61] Let f : S → T be a measurable mapping, let η be a measurable measure on S, and let g be a measurable function on T . Then ´ T gdf#η = ´ S (g ◦ f)dη.
Lemma 3 [Lemma 6 in Paty & Cuturi (2019)] For all φ, ψ and µ ∈ P(X ), ν ∈ P(Y), one has Π(φ#µ, ψ#ν) = {(φ ⊗ ψ)#γ s.t. γ ∈ Π(µ, ν)}, where φ ⊗ ψ : X × Y → X × Y such that (φ⊗ ψ(x, y) = (φ(x), ψ(y)) for all x, y ∈ X × Y.
• (i)HWDr(µ, µ) is finite. In one hand, we assume that µ ∈Pr(X ) and ν ∈Pr(Y), hence its r-th moments are finite, i.e., Mr(µ) = ( ´ X ‖x‖ rdµ(x) )1/r <∞ and Mr(ν) = ( ´ Y ‖y‖ rdν(y) )1/r < ∞. In the other hand, the following holds for all parameter θ ∈ Sd−1 and a couple (φ, ψ)-embeddings,
Wrr (µφ,θ, νψ,θ) = inf π∈Π(Pφ(θ)#µ,ψ(θ)#ν) ˆ R×R |u− u′|rdπ(u, u′)
(Lemma 3) = inf
γ∈Π(µ,ν) ˆ X×Y |φ(θ)>x− ψ(θ)>y|rdγ(x, y)
≤ 2r−1 inf γ∈Π(µ,ν) ˆ X×Y ( |φ(θ)>x|r + |ψ(θ)>y|r ) dγ(x, y)
= 2r−1 inf γ∈Π(µ,ν) (ˆ X |φ(θ)>x|rdµ(x) + ˆ Y |ψ(θ)>y|rdν(y) ) ,
where we use the facts that (s+ t)r ≤ 2r−1(sr + tr),∀s, t ∈ R+ and that any γ transport plan has marginals µ on X and ν on Y . By Cauchy–Schwarz inequality, we get
Wrr (µφ,θ, νψ,θ) ≤ 2r−1 (ˆ X ‖φ(θ)‖r‖x‖rdµ(x) + ˆ Y ‖ψ(θ)‖r‖y‖rdν(y) ) = 2r−1 ( Mrr (µ) +M r r (ν) ) .
Then, ( Eθ∼Θ [ Wrr (µφ,θ, νψ,θ) ])1/r ≤ 2 r−1r (Mrr (µ) + Mrr (ν))1/r ≤ 2 r−1r (Mr(µ) + Mr(ν)). Finally, one has thatHWDr(µ, ν) ≤ 2 r−1 r ( Mr(µ) +Mr(ν) ) .
• (ii) Non-negativity and symmetry. Together the non-negativity, symmetry of Wasserstein distance and the decoupling property of iterated infima (or principle of the iterated infima) yield the nonnegativity and symmetry of the distributional sliced sub-embedding distance.
• (ii) HWDr(µ, µ) = 0. Let φ and φ′ two embeddings for projecting the same distribution µ. Without loss of generality, we suppose that the corresponding (φ, φ′)-admissible constants C ′φ ≤ Cφ, hence MCφ′ ⊆MCφ . Using the fact that sup(A ∩B) ≤ supA ∧ supB, (with a ∧ b) = min(a, b)), se have, straightforwardly,
HWDr(µ, µ) = inf φ,φ′ sup Θ∈MCφ∩MCφ′
( Eθ∼Θ [ Wrr (µφ,θ, µφ′,θ) ]) 1r ≤ inf φ,φ′ ( sup
Θ∈MCφ
( Eθ∼Θ [ Wrr (µφ,θ, µφ′,θ) ]) 1r ∧ sup Θ∈MC
φ′
( Eθ∼Θ [ Wrr (µφ,θ, µφ′,θ) ]) 1r)
= inf φ inf φ′
( sup
Θ∈MCφ
( Eθ∼Θ [ Wrr (µφ,θ, µφ′,θ) ]) 1r ∧ sup Θ∈MC
φ′
( Eθ∼Θ [ Wrr (µφ,θ, µφ′,θ) ]) 1r)
≤ inf φ
( sup
Θ∈MCφ
( Eθ∼Θ [ Wrr (µφ,θ, µφ,θ) ]) 1r ∧ sup Θ∈MC
φ′
( Eθ∼Θ [ Wrr (µφ,θ, µφ,θ) ]) 1r)
≤ inf φ sup Θ∈MCφ
( Eθ∼Θ [ Wrr (µφ,θ, µφ,θ) ]) 1r = 0.
• (iii) One has (
1 d ) 1 r infφ,ψ maxθ∈Sd−1Wr(µφ,θ, νψ,θ) ≤ HWDr(µ, ν) ≤ infφ,ψ maxθ∈Sd−1Wr(µφ,θ, νψ,θ).
Since MCφ ∩MCψ ⊂M1 andWrr (µφ,θ, νψ,θ) ≤ maxθ∈Sd−1Wrr (µφ,θ, νψ,θ) we find that
sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr (µφ,θ, νψ,θ) ]) 1r ≤ sup Θ∈M1 ( Eθ∼Θ [ Wrr (µφ,θ, νψ,θ) ]) 1r ≤ (
max θ∈Sd−1
Wrr (µφ,θ, νψ,θ) )1/r
≤ max θ∈Sd−1 Wr(µφ,θ, νψ,θ),
which entails that HWDr(µ, ν) ≤ infφ,ψ maxθ∈Sd−1Wr(µφ,θ, νψ,θ). Moreover, since the (φ, ψ)admissible constants Cφ and Cψ satisfy Cφ ≥ Uφd and Cψ ≥ Uψ d hence σ̄ = ∑d l=1 1 dδθl ∈MCφ ∩ MCψ , where we set θ1 = argmaxθ∈Sd−1Wr(µφ,θ, νψ,θ). We then obtain
HWDr(µ, ν) ≥ inf φ,ψ
( Eθ∼σ̄ [ Wrr (µφ,θ, νψ,θ) ]) 1r = inf φ,ψ ( d∑ l=1 1 d Wrr (µφ,θl , νψ,θl) ) 1 r
≥ (1 d )1/r inf φ,ψ Wr(µφ,θ1 , νψ,θ1)
= (1 d )1/r inf φ,ψ max θ∈Sd−1 Wr(µφ,θ, νψ,θ).
• (iv) For p = q, HWD is upper bound by the distributional Wasserstein distance (DSW) . Let us first recall the DSW distance: let C > 0 and set MC = {Θ ∈P(Sd−1) : Eθ,θ′∼Θ[|〈θ, θ′〉|] ≤ C}.
DSWr(µ, ν) = sup Θ∈MC
( Eθ∼Θ [ Wrr (µθ, νθ) ]) 1r .
We have that the case of a identity couple of embeddings, φ = Id, ψ = Id, the probability measure set MCφ ,MCψ = MC , then it is trivial thatHWDr(µ, ν) ≤ DSWr(µ, ν).
• (v) Rotation invariance. Note that (R#µ)φ,θ = Pφ(θ)#(R#µ) = (Pφ(θ) ◦ R)#µ, and for all x ∈ Rp, using the adjoint operator R∗, (R∗ = R−1), (Pφ(θ) ◦ R)(x) = 〈φ(θ), R(x)〉 = 〈R∗(φ(θ)), x〉 = PR∗◦φ(θ)(x). Then, (R#µ)φ,θ = (PR∗◦φ(θ))#µ. Analogously, one has (Q#ν)ψ,θ = (PQ∗◦ψ(θ))#ν. Moreover,
MCφ = { Θ ∈P(Sd−1) : Eθ,θ′∼Θ[|〈φ(θ), φ(θ′)〉|] }
= { Θ ∈P(Sd−1) : Eθ,θ′∼Θ[|〈(R∗ ◦ φ)(θ), (R∗ ◦ φ)(θ′)〉|] }
= MCR∗◦φ .
Then we have similarly MCψ = MCQ∗◦ψ . This implies
HWDr(R#µ,Q#ν) = inf φ,ψ sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr ((R#µ)φ,θ, (Q#ν)ψ,θ) ]) 1r = inf φ,ψ
sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr ((PR∗◦φ(θ))#µ, PQ∗◦ψ(θ))#ν) ]) 1r = inf φ,ψ
sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr (µR∗◦φ,θ, µR∗◦φ,θ ]) 1r = inf φ,ψ
sup Θ∈MCR∗◦φ∩MCQ∗◦ψ
( Eθ∼Θ [ Wrr (µR∗◦φ,θ, νQ∗◦φ,θ ]) 1r = inf φ,ψ
sup Θ∈MCR∗◦φ∩MCQ∗◦ψ
( Eθ∼Θ [ Wrr (µR∗◦φ,θ, νQ∗◦ψ,θ ]) 1 r
= inf φ′=R∗◦φ,ψ′=Q∗◦ψ sup Θ∈MC
φ′ ∩MC ψ′
( Eθ∼Θ [ Wrr (µφ′,θ, νψ′,θ ]) 1r = HWDr(µ, ν).
• (vi) Translation quasi-invariance. We have
HWDr(Tα#µ, Tβ#ν) = inf φ,ψ sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr ((Tα#µ)φ,θ, (Tβ#ν)ψ,θ) ]) 1r .
By Lemmas 3 and 2 , we have Wrr ((Tα#µ)φ,θ, (Tβ#ν)ψ,θ
= inf γ∈Π((Tα#µ)φ,θ,(Tβ#ν)ψ,θ)) ˆ R2 |u− v|rdγ(u, v)
= inf γ∈Π((Pφ(θ)◦Tα)#µ,(Pψ(θ)◦Tβ)#ν) ˆ R2 |u− v|rdγ(u, v)
= inf γ∈Π(µ,ν) ˆ X×Y |Pφ(θ) ◦ Tα(x)− Pψ(θ) ◦ Tβ)(y)|rdγ(x, y)
= inf γ∈Π(µ,ν) ˆ X×Y |(Pφ(θ)(x)− Pψ(θ)(y)) + (Pφ(θ)(α)− Pψ(θ)(β))|rdγ(x, y)
≤ 2r−1 (
inf γ∈Π(µ,ν) ˆ X×Y |(Pφ(θ)(x)− Pψ(θ)(y))|rdγ(x, y) + |Pφ(θ)(α)− Pψ(θ)(β))|r )
≤ 2r−1 (
inf γ∈Π(µ,ν) ˆ X×Y |(Pφ(θ)(x)− Pψ(θ)(y))|rdγ(x, y) + (‖α‖+ ‖β‖)r ) .
Thanks to Minkowski inequality,
sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr ((Tα#µ) φ θ , (Tβ#ν) ψ θ ) ]) 1r
≤ 2r−1 sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ inf
γ∈Π(µ,ν) ˆ X×Y |(Pφ(θ)(x)− Pψ(θ)(y))|rdγ(x, y) ]) 1 r
+ 2r−1(‖α‖+ ‖β‖) sup Θ∈MCφ∩MCψ
( Θ(Sd−1) ) 1 r
≤ 2r−1 sup Θ∈MCφ∩MCψ
( Eθ∼Θ [ Wrr (µφ,θ, νψ,θ) ]) 1r + 2r−1(‖α‖+ ‖β‖).
Therefore, we getHWDr(Tα#µ, Tβ#ν) ≤ 2r−1HWDr(µ, ν) + 2r−1(‖α‖+ ‖β‖).
B IMPLEMENTATION
This section graphically describes the learning procedure in Algorithm 1. It also provides the training details not exposed in the main body of the paper.
B.1 LEARNING SCHEME
We present in Figure 7 the updated graphics of our approach, highlighting the main components : the distributional part is ensured by a first deep neural network as is each of the mappings. As each of the networks should be learned, we included the part of the loss functions associated with each network (blue fonts correspond to minimization, whereas red fonts correspond to maximization, see Algorithm section).
B.2 TRAINING DETAILS
Our experimental evaluations on shape datasets for scalability contrast GW, SGW and HWD. Regarding or classification under isometry transformations, we additionnally consider RI-SGW. Used
hyper-parameters for those experiments are detailed below. Notice that SGW, RI-SGW and HWD rely on K, the number of projections sampled uniformly over the unit sphere. This K may vary from a method to another.
1. SGW: K. 2. RI-SGW: λRI-SGW, the learning rate and T , the maximal number of iterations for solving 4
over the Stiefeld manifold. 3. HWD: beyond K and the latent space dimension d, it requires the parametrization of φ, ψ
and f as deep neural networks and their optimizers. For solving the min-max problem by an alternating optimization scheme we use N inner loops and T number of epochs.
For SGW and RI-SGW we use the code made available by their authors and cite the related reference Vayer et al. (2018) as they require. We use POT toolbox Flamary et al. (2021) to compute GW distance.
Scalability This experiment measures the average running time to compute OT-based distance between two pairs of shapes made of n 3D-vertices. 100 pairs of shapes were considered and n varies in {100, 250, 500, 1000, 1500, 2000}. We choose KSGW = 1000 (as a default value).
For HWD, the mapping function f is designed as a deep network with 2 dense hidden layers of size 50. Regarding both φ and ψ, they have also the same architecture as f (with adapted input and output layers) but the hidden layers are 10-dimensional. Adam optimizers with default parameters are selected to train them. Finally we consider KHWD = 10, d = 5, T = 50, N = 1 as default values. Notice also that the regularization parameters λC and λa are set to 1.
The used ground cost distance for GW distance is the geodesic distance.
Classification under transformations invariance For this experiment, we consider the same set of hyper-parameters as for Scalability evaluation on shape datasets. Besides, the supplementary competitor RI-SGW was trained by setting KRI-SGW = 1000 = KSGW = 1000, λRI-SGW = 0.01, TRI-SGW = 500. Notice that due to the high-resolution of the meshes (more than 19K threedimensional vertices), RI-SGW and GW were not able to produce the pairwise-distance matrix used in 1NN classification after several hours.

heading: C ADDITIONAL EXPERIMENTAL RESULTS

Review summaries:
This paper introduces the so-called Heterogeneous Wasserstein Discrepancy (HWD) between two probability measures supported on Euclidean spaces of different dimensions, a case scenario appearing in several applications. The HWD is built on ideas from the sliced Wasserstein distance and measure embedding : the two distributions are projected into a one-dimensional space, in which the classical Wasserstein distance (closed-form in 1D) is computed. Moreover, an algorithm to compute such a divergence is presented, and its usefulness is illustrated on synthetic examples, validated on real-world data and compared to standard others methods.

The authors propose a variant of Wasserstein for distributions supported on different Euclidean spaces (e.g., different dimensions). The main ideas are to combine the embedding (for distributions supported in different spaces -- Alaya et al. 2020) and distributional sliced Wasserstein (Nguyen et al., 2020) to generate random slicing projections into 1d-space, and leverage the 1d-Wasserstein for the comparison. The authors apply the proposed method for generating modeling and in query framework. 

In the paper, the authors proposed heterogeneous Wasserstein discrepancy (HWD) to compare incomparable distributions. The main idea of HWD is via distributional slicing, which had been employed in the work of Nguyen et al. [2020].