The Unbalanced Gromov Wasserstein Distance: Conic Formulation and Relaxation
Year: 2022

Abstract:
Comparing metric measure spaces (i.e. a metric space endowed with a probability distribution) is at the heart of many machine learning problems. The most popular distance between such metric measure spaces is the Gromov-Wasserstein (GW) distance, which is the solution of a quadratic assignment problem. The GW distance is however limited to the comparison of metric measure spaces endowed with a probability distribution. To alleviate this issue, we introduce two Unbalanced Gromov-Wasserstein formulations: a distance and a more tractable upper-bounding relaxation. They both allow the comparison of metric spaces equipped with arbitrary positive measures up to isometries. The first formulation is a positive and definite divergence based on a relaxation of the mass conservation constraint using a novel type of quadratically-homogeneous divergence. This divergence works hand in hand with the entropic regularization approach which is popular to solve large scale optimal transport problems. We show that the underlying non-convex optimization problem can be efficiently tackled using a highly parallelizable and GPU-friendly iterative scheme. The second formulation is a distance between mm-spaces up to isometries based on a conic lifting. Lastly, we provide numerical experiments on synthetic examples and domain adaptation data with a Positive-Unlabeled learning task to highlight the salient features of the unbalanced divergence and its potential applications in ML.

heading: 1 Introduction
text: Comparing data distributions on different metric spaces is a basic problem in machine learning. This class of problems is for instance at the heart of surfaces [Bronstein et al., 2006] or graph matching [Xu et al., 2019] (equipping the surface or graph with its associated geodesic distance), regression problems in quantum chemistry [Gilmer et al., 2017] (viewing the molecules as distributions of points in R3) and natural language processing [Grave et al., 2019, Alvarez-Melis and Jaakkola, 2018] (where texts in different languages are embedded as points distributions in different vector spaces).
Metric measure spaces. The mathematical way to formalize these problems is to model the data as metric measure spaces (mm-spaces). A mm-space is denoted as X = (X, d, µ) where X is a complete separable set endowed with a distance d and a positive Borel measure µ ∈M+(X). For instance, if X = (xi)i is a finite set of points, then µ = ∑ imiδxi (here δxi is the Dirac mass at xi) is simply a set of positive weights mi = µ({xi}) ≥ 0 associated to each point xi, which accounts for its mass or importance. For instance, setting some mi to 0 is equivalent to removing the point xi. We
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
refer to Sturm [2012] for a mathematical account on the theory of mm-spaces. In all the applications highlighted above, it makes sense to perform the comparisons up to isometric transformations of the data. Two mm-spaces X = (X, dX , µ) and Y = (Y, dY , ν) are considered to be equal (denoted X ∼ Y) if they are isometric, meaning that there is a bijection ψ : spt(µ)→ spt(ν) (where spt(µ) is the support of µ) such that dX(x, y) = dY (ψ(x), ψ(y)) and ψ]µ = ν. Here ψ] is the push-forward operator, so that ψ]µ = ν is equivalent to imposing ν(A) = µ(ψ−1(A)) for any set A ⊂ Y . For discrete spaces where µ = ∑ imiδxi , then one should have ν = ψ]µ = ∑ imiδψ(xi). As highlighted by Mémoli [2011], considering mm-spaces up to isometry is a powerful way to formalize and analyze a wide variety of problems such as matching, regression and classification of distributions of points belonging to different spaces. Most often, the objects of interest come with a natural distance such as an intrinsic or extrinsic distance and the uniform measure is the usual choice to make mm-spaces widely applicable. The key to unlock all these problems is the computation of a distance between mm-spaces up to isometry. So far, existing distances (reviewed below) assume that µ is a probability distribution, i.e. µ(X) = 1. This constraint is not natural and sometimes problematic for most of the practical applications to machine learning. The goal of this paper is to alleviate this restriction. We define for the first time a class of distances between unbalanced metric measure spaces, these distances being upper-bounded by divergences which can be approximated by an efficient numerical scheme.
Csiszár divergences The simplest case is when X = Y and one simply ignores the underlying metric. One can then use Csiszár divergences (or ϕ-divergences), which perform a pointwise comparison (in contrast with optimal transport distances, which perform a displacement comparison). It is defined using an entropy function ϕ : R+ → [0,+∞], which is a convex, lower semi-continuous, positive function with ϕ(1) = 0. The Csiszár ϕ-divergence reads Dϕ(µ|ν) , ∫ X ϕ ( dµ dν ) dν +
ϕ′∞ ∫ X dµ⊥, where µ = dµdν ν + µ ⊥ is called the Radon-Nikodym or the Lebesgue decomposition of µ with respect to ν and ϕ′∞ = limr→∞ ϕ(r)/r ∈ R ∪ {+∞} is called the recession constant. This divergence Dϕ is convex, positive, 1-homogeneous and weak* lower-semicontinuous, see Liero et al. [2015] for details. Particular instances of ϕ-divergences are Kullback-Leibler (KL) for ϕ(r) = r log(r)− r + 1 (note that ϕ′∞ =∞) and Total Variation (TV) for ϕ(r) = |r − 1|.
Balanced and unbalanced optimal transport. If the common embedding space X is equipped with a distance d(x, y), one can use more elaborated methods such as optimal transport (OT) distances, which are computed by solving convex optimization problems. This type of methods has proven useful for ML problems as diverse as domain adaptation [Courty et al., 2014], supervised learning over histograms [Frogner et al., 2015] and unsupervised learning of generative models [Arjovsky et al., 2017]. In this case, the extension from probability distributions to arbitrary positive measures (µ, ν) ∈ M+(X)2 is now well understood and corresponds to the theory of unbalanced OT. Following Liero et al. [2015], Chizat et al. [2018a], a family of unbalanced Wasserstein distances is defined by solving
UW(µ, ν)q , inf π∈M(X×X)
∫ λ(d(x, y))dπ(x, y) + Dϕ(π1|µ) + Dϕ(π2|µ). (1)
Here (π1, π2) are the two marginals of the joint distribution π, defined by π1(A) = π(A× Y ) for A ⊂ X . The mapping λ : R+ → R and exponent q ≥ 1 should be chosen wisely to ensure for instance that UW defines a distance (see Section 2.2.1). It is frequent to take ρDϕ instead of Dϕ (i.e. take ψ = ρϕ) to adjust the strength of the marginals’ penalization. Balanced OT is retrieved with the convex indicator ϕ = ι{1} (i.e. ϕ(1) = 0 and ϕ(x) = +∞ otherwise) or by taking the limit ρ → +∞, which enforces π1 = µ and π2 = ν. When 0 < ρ < +∞, unbalanced OT operates a trade-off between transportation and creation of mass, which is crucial to be robust to outliers in the data and to cope with mass variations in the modes of the distributions. For supervised tasks, the value of ρ should be cross-validated to obtain the best performances. Its use is gaining popularity in applications, such as medical imaging registration [Feydy et al., 2019a], videos [Lee et al., 2019], generative learning [Balaji et al., 2020] and gradient flow to train neural networks [Chizat and Bach, 2018, Rotskoff et al., 2019]. Furthermore, existing efficient algorithms for balanced OT extend to this unbalanced problem. In particular Sinkhorn’s iterations, introduced in ML for balanced OT by Cuturi [2013], extend to unbalanced OT [Chizat et al., 2018b, Séjourné et al., 2019], as detailed in Section 3.
The Gromov-Wasserstein distance and its applications. The Gromov-Wasserstein (GW) distance [Mémoli, 2011, Sturm, 2012] generalizes the notion of OT to the setting of mm-spaces up to
isometries. It replaces the linear cost ∫ λ(d)dπ of OT by a quadratic function. It reads
GW(X ,Y)q , min π∈M+(X×Y )
{∫ λ(|dX(x, x′)− dY (y, y′)|)dπ(x, y)dπ(x′, y′) :
π1 = µ π2 = ν
} . (2)
It is proved in Mémoli [2011], Sturm [2012] that GW defines with λ(t) = tq a distance up to isometries on balanced mm-spaces (i.e. the measures are probability distributions). The GW distance is applied successfully in natural language processing for unsupervised translation learning [Grave et al., 2019, Alvarez-Melis and Jaakkola, 2018], in generative learning for objects lying in spaces of different dimensions [Bunne et al., 2019] and to build VAE for graphs [Xu et al., 2020]. It has been adapted for domain adaptation over different spaces [Redko et al., 2020]. It is also a relevant distance to compute barycenters between graphs or shapes [Vayer et al., 2018, Chowdhury and Needham, 2020]. When (X ,Y) are Euclidean spaces, this distance compares distributions up to rigid isometry, and is closely related (but not equal) to metrics defined by procrustes analysis [Grave et al., 2019, Alvarez-Melis et al., 2019]. The problem (2) is non convex because the quadratic form∫ λ(|dX − dY |)dπ ⊗ π is not positive in general. It is in fact closely related to quadratic assignment problems [Burkard et al., 1998], which are used for graph matching problems, and are known to be NP-hard in general. Nevertheless, non-convex optimization methods have been shown to be successful in practice to use GW distances for ML problems. This includes for instance alternating minimization [Mémoli, 2011, Redko et al., 2020] and entropic regularization [Peyré et al., 2016, Gold and Rangarajan, 1996].
Related works and contributions. The concomitant work of De Ponti and Mondino [2020] extends the Lp transportation distance defined in Sturm et al. [2006] to unbalanced mm-spaces and studies its geometric properties. This distortion distance is not equivalent to the GW distance, and is more difficult to estimate numerically because it explicitly imposes a triangle inequality constraint in the optimization problem. The work of Chapel et al. [2020] relaxes the GW distance to the unbalanced setting by hybridizing GW with partial OT [Figalli, 2010] for unsupervised labeling. It ressembles one particular setting of our formulation, but with some important differences, detailed in Section 2. Our construction is also connected to partial matching methods, which find numerous applications in graphics and vision [Cosmo et al., 2016]. In particular, Rodola et al. [2012] introduces a mass conservation relaxation of the GW problem.
The two main contributions of this paper are the definition of two formulations relaxing the GW distance. The first one is called the Unbalanced Gromov-Wasserstein (UGW) divergence and can be computed efficiently on GPUs. The second one is called the Conic Gromov-Wasserstein distance (CGW). It is proved to be a distance between mm-spaces endowed with positive measures up to isometries, as stated in Theorem 1 which is the main theoretical result of this paper. We also prove in Theorem 1 that UGW can be used as a surrogate upper-bounding CGW. We present those concepts and their properties in Section 2. We also detail in Section 3 an efficient computational scheme for a particular setting of UGW. This method computes an approximate stationary point of a biconvex relaxation of our formulations. Even though it is a lower bound of the original problem, we provide in Theorem 3 conditions ensuring the tightness of this relaxation in many cases of interest. The algorithm leverages the strength of entropic regularization and the Sinkhorn algorithm, namely that it is GPU-friendly and defines smooth loss functions amenable to back-propagation for ML applications. Section 4 provides some numerical experiments to highlight the qualitative behavior of this algorithm and its ability to cope with outliers and mass variations in the modes of the distributions. We illustrate numerically the tightness of the relation between UGW and CGW, showing that UGW is a reasonnable proxy of a distance, at least locally. We provide an application of our divergence in the positive unlabeled learning setting, using domain adaptation data, and display results which are at par or outperform the computable competitor Chapel et al. [2020].

heading: 2 Unbalanced Gromov-Wasserstein formulations
text: We present in this section our two new formulations and their properties. The first one, called UGW, is exploited in Sections 3 and 4 to derive an efficient algorithm used in numerical experiments. The second one, called CGW, defines a distance between mm-spaces up to isometries. Those results build upon the work of Liero et al. [2015], and a summary of the construction of UOT is detailed in Appendix A. In all what follows, we consider complete separable mm-spaces endowed with a metric and a positive measure.

heading: 2.1 The unbalanced Gromov-Wasserstein divergence
text: This new formulation makes use of quadratic ϕ-divergences, defined as D⊗ϕ (ρ|ν) , Dϕ(ρ⊗ ρ|ν⊗ ν), where ρ ⊗ ρ ∈ M+(X2) is the tensor product measure defined by d(ρ ⊗ ρ)(x, y) = dρ(x)dρ(y). Note that D⊗ϕ is not a convex function in general.
Definition 1 (Unbalanced GW). The Unbalanced Gromov-Wasserstein divergence is defined as UGW(X ,Y) = infπ∈M+(X×Y ) L(π) , G(π) + D⊗ϕ (π1|µ) + D⊗ϕ (π2|ν).
This definition can be understood as an hybridation between (1) and (2) but with a twist: one needs to use the quadratic divergence D⊗ϕ in place of Dϕ. To the best of our knowledge, it is the first time such quadratic divergences are being used and studied. In the TV case, this is the most important distinction between UGW and partial-GW [Chapel et al., 2020]. Note also that the balanced GW distance (2) is recovered as a particular case when using ϕ = ι{1} or by letting ρ → +∞ for an entropy ψ = ρϕ.
Using quadratic divergences results in UGW being 2-homogeneous: for θ ≥ 0, writing (Xθ,Yθ) equiped with (θµ, θν), one has θ−2UGW(Xθ,Yθ) = UGW(X ,Y). When using non tensorized ϕ-divergences, the resulting unbalanced Gromov-Wassertein functional between Xθ and Yθ have very different and inconsistent behaviors when θ → 0 and θ → +∞. Indeed, once normalized by θ−2 and θ−1, one obtains respectively balanced GW and a Hellinger-type distance. Using tensorized divergences ensures that the behavior does not depends on θ. It is also fundamental to connect UGW with our distance CGW, see Theorem 1 and Appendix B.
We first prove the existence of optimal plans π minimizing L, which holds for the three key settings of Section 2.2.1, namely for KL, TV, and for compact metric spaces (such as finite pointclouds and graphs). All proofs are deferred in Appendix B.
Proposition 1 (Existence of minimizers). We assume that (X,Y ) are compact and that either (i) ϕ superlinear, i.e ϕ′∞ = ∞, or (ii) λ has compact sublevel sets in R+ and 2ϕ′∞ + inf λ > 0. Then there exists π ∈M+(X × Y ) such that UGW(X ,Y) = L(π).
The following proposition ensures that the functional UGW can be used to compare mm-spaces.
Proposition 2 (Definiteness of UGW). Assume that ϕ−1({0}) = {1} and λ−1({0}) = {0}. Then UGW(X ,Y) ≥ 0 and is 0 if and only if X ∼ Y .
We end this section with a reformulation of UGW which is important to make the connection with the second formulation CGW of the following section. It splits UGW into two parts: the term ϕ(0)(|(µ⊗µ)⊥|+ |(ν⊗ν)⊥|) accounts for the pure creation/destruction of mass and a new transport cost Lc accounts for the remaining part (partial/pure transport and partial creation/destruction of mass).
Lemma 1. Defining Lc(a, b) , c + aϕ(1/a) + bϕ(1/b), and writing (f , dµdπ1 , g , dν dπ2 ) the Lebesgue densities of (µ, ν) w.r.t. (π1, π2) such that µ = fπ1 + µ⊥ and ν = gπ2 + ν⊥, one has
L(π) = ∫ X2×Y 2 Lλ(|dX−dY |)(f ⊗ f, g ⊗ g)dπdπ + ϕ(0)(|(µ⊗ µ) ⊥|+ |(ν ⊗ ν)⊥|). (3)
Proof. Write f = dµdπ1 and g = dν dπ2 . The Lebesgue decompositions read µ ⊗ µ = (f ⊗ f)π1 ⊗ π1 + (µ ⊗ µ)⊥ and ν ⊗ ν = (g ⊗ g)π2 ⊗ π2 + (ν ⊗ ν)⊥, thanks to the tensorized structure of the decomposed plans. To prove Equation (3), we need to define the reverse entropy Liero et al. [2015] such that Dϕ(α|µ) = Dψ(µ|α), where ψ(x) , xϕ( 1x ) is also an entropy function satisfying ψ′∞ = ϕ(0). One then has
L(π) = ∫ X2×Y 2 λ(Γ)dπdπ + D⊗ϕ (π1|µ) + D⊗ϕ (π2|ν)
= ∫ X2×Y 2 λ(Γ)dπdπ + D⊗ψ (µ|π1) + D ⊗ ψ (ν|π2)
L(π) = ∫ X2×Y 2 λ(Γ)dπdπ + ∫ X2 ψ(f ⊗ f)dπ1dπ1 + ∫ Y 2 ψ(g ⊗ g)dπ2dπ2
+ ϕ(0)(|(µ⊗ µ)⊥|+ |(ν ⊗ ν)⊥|)
= ∫ X2×Y 2 Lλ(Γ)(f ⊗ f, g ⊗ g)dπdπ + ϕ(0)(|(µ⊗ µ)⊥|+ |(ν ⊗ ν)⊥|).
Using the definition of ψ in Lc ends the proof.

heading: 2.2 The conic Gromov-Wasserstein distance
text: We introduce a second “conic” formulation of unbalanced GW, which is connected to UGW, and whose construction is inspired by the conic formulation of UOT (see Appendix A for an overview).

heading: 2.2.1 Background on cone sets and distances
text: The conic formulation lifts a point x ∈ X to a couple (x, r) ∈ X × R+ where r encodes some (power of a) mass. Then we seek optimal transport plans defined over C[X] , X × R+/(X × {0}), where coordinates (x, r = 0) with no mass are merged into a single point 0X called the apex of the cone. In the sequel, points of X × R+ are noted (x, r), while [x, r] are quotiented points of C[X]. While transport plans depend on variables ([x, r], [y, s]) and ([x′, r′], [y′, s′]) in C[X]×C[Y ], the transportation cost involved in our conic formulation only makes use of the 2-D cone C[R+] over R+ endowed with the distance |u−v| (note that any other distance on R could be used as well). More specifically, we consider coordinates of the form ([u, a], [v, b]) = ([dX(x, x′), rr′], [dY (y, y′), ss′]) ∈ C[R+] × C[R+]. Thus we now describe conic discrepancies D on C[R+], which are defined for (p, q) ≥ 1 as D([u, a], [v, b])q , Hλ(|u−v|)(ap, bp), where Hc(ap, bp) , infθ≥0 θLc(a p θ , bp
θ ) is the perspective transform of Lc introduced in Lemma 1. The intuition underpinning the definition of this cost is that the perspective transform accounts for the possibility to rescale a transport plan π by a scalar θ but the scaling is performed pointwise instead of globally. In general D is not a distance, but it is always definite as stated by this result proved in Appendix A. Proposition 3. Assume λ−1({0}) = {0}, ϕ−1({0}) = {1} and ϕ is coercive. Then D is definite on C[R+], i.e. D([u, a], [v, b]) = 0 if and only if (a = b = 0) or (a = b and u = v).
Of particular interest are those ϕ where D is a distance, which necessitates a careful choice of λ, p and q. We now detail three examples where this is the case. Gaussian Hellinger distance (GH). When Dϕ = KL, λ(t) = t2 and q = p = 2, then one has D([u, a], [v, b])2 = a2 + b2 − 2abe−|u−v|/2. This cone distance [Burago et al., 2001] is further generalized by De Ponti [2019] who shows that D is a distance for power entropies ϕ(s) = s
p−p(s−1)−1 p(p−1)
if p ≥ 1 (the case p = 1 corresponding to Dϕ = KL). Hellinger-Kantorovich (HK) / Wasserstein-Fisher-Rao distance (WFR). When Dϕ = KL, λ(t) = − log cos2(t∧ π2 ) and q = p = 2, then one hasD([u, a], [v, b])
2 = a2 +b2−2ab cos(π2 ∧|u− v|). This construction, which might seem peculiar, corresponds to the one used to make unbalanced OT a geodesic distance, as detailed in [Liero et al., 2015, Chizat et al., 2018a]. Partial optimal transport distance (PT). When Dϕ = TV, λ(t) = tq, q ≥ 1 and p = 1, then D([u, a], [v, b])q = a+ b− (a ∧ b)(2− |u− v|q)+ defines a cone distance [Chizat et al., 2018a].

heading: 2.2.2 Definitions and properties
text: The conic formulation consists in solving a GW problem on the cone, with the addition of two linear constraints. Informally speaking, Lc from Lemma 1 becomes D, the term (|(µ⊗ µ)⊥|+ |(ν ⊗ ν)⊥|) is taken into account by the constraints (5) below, and the variables (f, g) are replaced by (rp, sp). It reads CGW(X ,Y) , infα∈Up(µ,ν)H(α) where
H(α) , ∫ D([dX(x, x′), rr′], [dY (y, y′), ss′])q dα([x, r], [y, s])dα([x′, r′], [y′, s′]), (4)
and Up(µ, ν) is defined as the set Up(µ, ν) , { α ∈M+(C[X]× C[Y ]), ∫ R+ rpdα1(·, r) = µ, ∫ R+ spdα2(·, s) = ν } . (5)
It is similar to the conic formulation of UW, see Appendix A. Note that similarly to the GW formulation (2) – and in sharp contrast with the conic formulation of UW – here the transport plans are defined on the cone C[X]× C[Y ] but the cost D is a distance on C[R+]. We present now the main contributions of this paper, proved in Appendix C. We state that CGW defines a distance under conditions that hold for the settings of Section 2.2.1, and that it is upperbounded by UGW. The divergence UGW can be approximated with efficient numerical schemes as detailed in Section 3.
Theorem 1. (i) The divergence CGW is symmetric, positive and definite up to isometries. (ii) If D is a distance on C[R+], then CGW1/q is a distance on the set of mm-spaces up to isometries. (iii) For any (Dϕ, λ, p, q) with associated cost D on the cone, one has UGW ≥ CGW.

heading: 3 Algorithms
text: We focus in this section on the numerical computation of the upper bound UGW using a bi-convex relaxation and derive an alternate minimization scheme coupled with entropic regularization. We also propose to approximate CGW by doing a similar alternate minimization, as detailed in Appendix E. We provide guarantees of tightness on the bi-convex relaxation for CGW (see Theorem 3). The computation of the distance CGW is heavy in practice because it requires an optimization over a lifted conic space, which needs to be discretized. Thus it does not scale to large problem for CGW, but allows to explore numerically how tight is the upper bound UGW ≥ CGW, see Section 4. The algorithm for UGW is presented on arbitrary measures, the special case of discrete measures being a particular case. The discretized formulas and algorithms are detailed in Appendix D, see also Chizat et al. [2018b], Peyré et al. [2016]. All implementations are available at https: //github.com/thibsej/unbalanced_gromov_wasserstein, and installable in Python with the command pip install unbalancedgw.

heading: 3.1 Bi-convex relaxation and tightness
text: In order to derive a simple numerical approximation scheme, following Mémoli [2011], we introduce a lower bound obtained by introducing two transportation plans. To further accelerate the method and enable GPU-friendly iterations, similarly to Gold et al. [1996], Solomon et al. [2016], we consider an entropic regularization. It reads, for any ε ≥ 0,
UGWε(X ,Y) , inf π L(π) + εKL⊗(π|µ⊗ ν) ≥ inf π,γ F(π, γ) + εKL(π ⊗ γ|(µ⊗ ν)⊗2), (6) and F(π, γ) , ∫ X2×Y 2 λ(|dX − dY |)dπ ⊗ γ + Dϕ(π1 ⊗ γ1|µ⊗ µ) + Dϕ(π2 ⊗ γ2|ν ⊗ ν) ,
where (γ1, γ2) denote the marginals of the plan γ. In the sequel we write Fε = F + εKL⊗. Note that in contrast to the entropic regularization of GW Peyré et al. [2016], here we use a tensorized entropy to maintain the overall homogeneity of the energy. A simple method to approximate this lower bound is to perform an alternate minimization on π and γ, which is known to converge for smooth ϕ to a stationary point since the coupling term in the functional is smooth [Tseng, 2001]. Note that if π ⊗ γ is optimal then so is (sπ)⊗ ( 1sγ) with s ≥ 0. Thus without loss of generality we can optimize under the constraint m(π) = m(γ) by setting s = √ m(γ)/m(π).
We now discuss the tightness of the bi-convex relaxation by generalizing a result of Konno. We first present a result which applies to general quadratic assignment problems, then state its application to our setting.
Theorem 2 (Tight relaxation). Let B a Banach space, let f : B 7→ R ∪ {+∞} be a function and let L : C ⊂ B 7→ R the function defined on the convex set C ⊂ B by L(π) = 12 〈π, k(π)〉+2f(π) where k is a symmetric bilinear map which is negative (not necessarily definite) on ∆C , Span({π − γ ; (π, γ) ∈ C}), that is, for any z ∈ ∆C, 〈z, kz〉 ≤ 0. Assume that there exists π0 ∈ C such that L(π0) < +∞, and define F(π, γ) , 12 〈π, k(γ)〉 + f(π) + f(γ). Then, for any (π∗, γ∗) ∈ arg minF(π, γ), we have F(π∗, π∗) = F(γ∗, γ∗) = F(π∗, γ∗). Moreover, if one assumes either that k is a definite kernel or f is strictly convex, one gets π∗ = γ∗.
The above Theorem 2 is proved in Appendix D. As an application, we now state our tightness result for GWε and CGW. In those settings the optimizers of the bi-convex relaxation are also optimal for the original problem. Theorem 3. For GWε with ε ≥ 0 or for CGW, assume that λ(t) = t2 and that (dX , dY ) are both conditionnally negative (or conditionally positive) kernels. Then the bi-convex relaxation of both problems is tight.
Proof. The proof for CGW is detailed in Appendix E, we prove the tightness for GWε. When λ(t) = t2 the kernel k = λ(|dX − dY |) is conditionally negative on the set {(π, γ), π1 = γ1 and π2 = γ2}, i.e. we have 〈(π−γ), k(π−γ)〉 ≤ 0 (see Maron and Lipman [2018]). For GWε one has π1 = γ1 = µ and π2 = γ2 = ν thanks to the constraints on marginals. Thus the kernel is negative semi-definite and the proof of Theorem 2 applies, hence the tightness of the relaxation.
Konno’s result Konno [1976] applies for unregularized (ε = 0), Balanced-GW. The novelty of Theorem 3 is its extension to both GWε and CGW. So far it is an open question whether the relaxation is tight or not for UGWε, because the above proof no longer holds. Note that in all our numerical simulations, our solvers always found solutions of UGWε such that π = γ when |dX − dY |2 is conditonally negative. The property that the kernel |dX − dY |2 is negative does not hold in general (e.g. for graph geodesic distances) and the tightness of the relaxation remains open in this setting. We know from [Maron and Lipman, 2018, Theorem 1] that it is conditionally negative semi-definite when both (dX , dY ) are conditionally negative kernels. Examples of distances which are negative kernels are tree metrics in the case of graphs, as well as Euclidean, spherical and hyperbolic distances over their respective manifolds Feragen et al. [2015]. In practice, when ε is small, we observed in the indefinite setting that the relaxation outputs more frequently spurious minima than in the negative semi-definite setting.

heading: 3.2 Alternate Sinkhorn minimization
text: Minimizing the lower bound (6) with respect to either π or γ is non-trivial for an arbitrary ϕ. We restrict our attention to the Kullback-Leibler case Dϕ = ρKL with ρ > 0, which can be addressed by solving a regularized and convex unbalanced problem as studied in Chizat et al. [2018b], Séjourné et al. [2019]. It is explained in the following proposition. Proposition 4. For a fixed γ, the optimal π ∈ arg min
π F(π, γ) + εKL(π ⊗ γ|(µ⊗ ν)⊗2) solves
min π
∫ cεγ(x, y)dπ(x, y) + ρm(γ)KL(π1|µ) + ρm(γ)KL(π2|ν) + εm(γ)KL(π|µ⊗ ν),
where m(γ) , γ(X × Y ) is the mass of γ, and where we define the cost associated to γ as cεγ(x, y), ∫ λ(|dX(x, ·)− dY (y, ·)|)dγ+ρ ∫ log(
dγ1 dµ )dγ1+ρ
∫ log(
dγ2 dν )dγ2+ε
∫ log( dγ
dµdν )dγ.
Algorithm 1 – UGW(X , Y , ρ, ε) Input: mm-spaces (X ,Y), relax. ρ, regul. ε Output: π, γ solving (6)
Init. π = γ = µ⊗ ν/ √ m(µ)m(ν), g = 0. while (π, γ) has not converged do Update π ← γ, then c← cεπ , ρ̃← m(π)ρ, ε̃← m(π)ε while (f, g) has not converged do f ← − ε̃ρ̃
ε̃+ρ̃ log
∫ e(g(y)−c(·,y))/ε̃dν(y)
g ← − ε̃ρ̃ ε̃+ρ̃
log ∫ e(f(x)−c(x,·))/ε̃dµ(x)
end while Upd. γ(x, y)←e f(x)+g(y)−c(x,y) ε̃ µ(x)ν(y)
Rescale γ ← √ m(π)/m(γ)γ
end while Return (π, γ).
Computing the cost cεγ for spaces X and Y of n points has in general a costO(n4) in time and memory. However, as explained for instance in Peyré et al. [2016], for the special case λ(t) = t2, this cost is reduced to O(n3) in time and O(n2) in memory. This is the setting we consider in the numerical simulations. This makes the method applicable for scales of the order of 104 points. For larger datasets one should use approximation schemes such as hierarchical approaches [Xu et al., 2019] or Nyström compression of the kernel [Altschuler et al., 2018].
The resulting alternate minimization method is detailed in Algorithm 1, see Appendix D for a discretized version. It uses the unbalanced Sinkhorn algorithm of Chizat et al. [2018b], Séjourné et al. [2019] as subiterations and takes π = µ ⊗ ν/ √ m(µ)m(ν) to initialize the updates. This Sinkhorn algorithm operates over a pair of continuous functions (so-called
Kantorovitch potentials) f(x) and g(y). For discrete spaces X and Y of size n, these functions are stored in vectors of size n, and that integral involved in the updates becomes a sum. Each iteration of Sinkhorn thus has a cost n2, and all the involved operation can be efficiently mapped to parallelizable GPU routines as detailed in Chizat et al. [2018b], Séjourné et al. [2019]. Another advantage of using an unbalanced Sinkhorn algorithm is its complexity O(n2/ε) to compute an ε-approximation, as stated in Pham et al. [2020], which should be compared to O(n2/ε2) operations for balanced Sinkhorn.
Note also that balanced GW is recovered as a special case when setting ρ → +∞, so that ρ̃/(ε̃ + ρ̃) → 1 should be used in the iterations. In order to speed up Sinkhorn inner-loops, especially for small values of ε, one can use linear extrapolation [Thibault et al., 2017] or non-linear Anderson acceleration [Anderson, 1965, Scieur et al., 2016].
There is an extra scaling step after computing γ involving the mass m(π). It corresponds to the scaling s of π⊗ γ such that m(π) = m(γ), and we observe that this scaling is key not only to impose this mass equality but also to stabilize the algorithm. Otherwise we observed that m(γ) < 1 < m(π) and underflows whenever m(γ)→ 0 and m(π)→∞.

heading: 4 Numerical experiments
text: This section presents simulations on synthetic examples to highlight the qualitative behavior of UGW and the tightness of the bound UGW ≥ CGW. Other illustrations on UGW are available in Appendix E. We end the section with a learning application of UGW in a positive-unlabeled setting, using domain adaptation data so as to compare with PGW Chapel et al. [2020]. In the synthetic experiments, µ and ν are probability distributions, which allows us to compare GW with UGW.
Robustness to imbalanced classes. In this first example, we take X = R3, Y = R2 and consider E2, E3, C and S to be uniform distributions on a 2D and 3D ellipse, a square and a sphere. We consider mm-spaces of different dimensions to emphasize the ability of (U)GW to compare different spaces. Figure 1 contrasts the transportation plan obtained by GW and UGW for a fixed µ = 0.5E3 + 0.5S and ν obtained using two different mixtures of E2 and C. The black segments show the largest entries of the transportation matrix π, for a sub-sampled set of points (to ease visibility), thus effectively displaying the matching induced by the plan. Furthermore, the width of the dots are scaled according to the mass of the marginals π1 ≈ µ and π2 ≈ ν, i.e. the smaller the point, the smaller is the amount of transported mass. This figure shows that the exact conservation of mass imposed by GW leads to a poor geometrical matching of the shapes which have different global mass. As this should be expected, UGW recovers coherent matchings. We suspect the alternate minimization algorithm is able to find the global minimum in these cases.
Tightness of the bound CGW≤UGW We propose to approximate CGW by doing a similar alternate minimization as for UGW, as detailed in Appendix E. This numerical scheme does not scale to large problems, but allows us to explore numerically how tight is the upper bound UGW ≥ CGW. Figure 2 highlights the fact that in Euclidean space X = Y = Rd, this bound seems to be tight when the two measures are sufficiently close. We consider discrete measures µ = 1n ∑ i δxi in
X = Y = Rd and νt = 1n ∑ i δyi where yi = xi + t∆i where ∆i are random perturbations and
denote (X ,Yt) the two mm-spaces associated to the Euclidean distance. As t → 0, µ and νt get closer, we observe numerically that UGW ≈ CGW. Figure 3 considers random points (xi)i and (yi)i and displays the histograms of the ratio CGW/UGW for n = 3. This shows that while the bound CGW ≤ UGW seems not tight, the ratio appears to be bounded even for points not being close. This numerical experiment suggests that UGW and CGW are locally equivalent and that UGW is in practice an acceptable proxy of the distance CGW. We leave for future works a tighter analysis of the gap between UGW and CGW.
their similar features, while negative samples would be ignored due to dissimilar features that induce a laziness to transport mass and match them.
We consider PU learning over the Caltech office dataset used for domain adaptation tasks (with domains Caltech (C) Griffin et al. [2007], Amazon (A), Webcam (W) and DSLR (D) Saenko et al. [2010]). The Caltech datasets are represented with two embeddings based on Surf and Decaf features Saenko et al. [2010], Donahue et al. [2014]. On the latter datasets, we perform PU learning over similar features (e.g. surf-C→ surf-* or decaf-C→ decaf-*) and from one feature format to the other (e.g. surf-C→ decaf-* or surf-C→ decaf-*). Those features are projected via PCA to subspaces of dimension 10 for surf features and 40 for decaf features. In the last task, one cannot use standard PU-method, and to the best of our knowledge, Unbalanced-GW methods are the only approaches for PU learning across different domains/features.
The procedure is the following. We solve the PU learning problem by computing the optimal plan π for UGW(X ,Y). We compute its first marginal π2 on Y , and predict the labels of some y ∈ Y as `(y) , sign(π2(y) − q) where q is the quantile of π2 corresponding to the proportion r of positives samples in Y . Following Chapel et al. [2020] which is adapted from Kato et al. [2018], Hsieh et al. [2019], this proportion r is assumed to be known. We report the accuracy of the prediction over the same 20 folds of the datasets, and use 20 other folds to validate the parameters of UGW. We consider 100 random samples for each fold of (X,Y ), a ratio of positive samples r = 0.1 for domains (C,A,W,D), and a ratio r = 0.2 for domains (C,A,W).
Since the GW objective is non-convex, the initialization of the minimization algorithms is key to obtain good perfor-
mances. In Chapel et al. [2020] and our experiments, for tasks where Y and Y belong to the same Euclidean space, (e.g. surf-C→ surf-*) we initialize π with the Partial-Wasserstein (PW) solution with a squared Euclidean cost. For cross-domain prediction (e.g. surf-C→ decaf-*), following Chapel et al. [2020], PGW is initialized with a list of plans built using a coarsened representation of the data with k-NN. While Chapel et al. [2020] makes use in an oracle manner of the plan providing the best accuracy, we modified their protocol and keep the plan which has the lowest PGW cost, which seems fairer, hence the difference in performance with Chapel et al. [2020]. To initialize UGW when X and Y do not belong to the same Euclidean space, we use a UOT solution of a matching between distance
histograms called FLB Mémoli [2011]. We define FLB in our UGW setting as FLB(X ,Y) , min ∫ X×Y |µ̄ ? dX − ν̄ ? dY |2dπ+ ρKL(π1|µ) + ρKL(π2|ν) + εKL(π|µ⊗ ν), (7)
where µ ? dX(x) , ∫ dX(x, x
′)dµ(x′) is the eccentricity, i.e. a histogram of aggregated distances, and µ̄ = µ/m(µ). Contrary to GW Mémoli [2011], there is a priori no link between FLB and UGW.
In the experiments we slightly generalize UGW and use two different marginal penalties ρ1KL
⊗(π1|µ) + ρ2KL⊗(π2|ν) with two parameters (ρ1, ρ2) to take into account shifts between domains/features. Note that PGW has a single parameter (which plays a role similar to (ρ1, ρ2)) which controls the cost of mass creation/destruction. We set ε = 2−9, which avoids introducing an extra parameter in the method. The value (ρ1, ρ2) ∈ {2−k, k ∈ J5, 10K}2 are cross validated for each task on the validation folds, and we report the average accuracy on the testing folds. We discuss in Appendix E the impact of reducing the number of parameters on the performance. Comparison with other methods – PU and PUSB Kato et al. [2018], Du Plessis et al. [2014] – are provided in Chapel et al. [2020] and we focus here on the comparison with PGW only.
The results are reported in Table 1. We display the performance of PGW, UGW and the initialization used for UGW to guarantee that using UGW does improve the performance. We observe that when the source and target dataset is the same (C→C tasks), the PW initialization performs better and PGW/UGW degrade the performance, so that in this setting Optimal Transport should be preferred over GW, which is to be expected. However when the domains are different, applying UGW improves the performance over the initialization (which is FLB) in almost all tasks. Note that in that case the methods PU, PUSB or PW cannot be used. Overall, this shows that GW methods are able to solve to some extent the PU learning problem across different spaces, and that using a “softer” KL penalties in UGW is at least at par with Partial GW, and performs better in some settings.

heading: 5 Conclusion and perspectives
text: This paper defines two Unbalanced Gromov-Wasserstein formulations: CGW and UGW. We prove that they are both positive and definite. We provide a scalable, GPU-friendly algorithm to compute UGW illustrate its applicability in learning tasks, and show that CGW is a distance between mmspaces up to isometry. These divergences and distances allow for the first time to blend in a seamless way the transportation geometry of GW with creation and destruction of mass. This hybridization is the key to unlock both theoretical and practical issues. This work opens new questions for future works, for instance removing the bias introduced by the use of entropic regularization, which is important for applications to ML. Note that such a debiasing was successfully applied for BalancedGW in Bunne et al. [2019] and is shown to lead to a valid divergence for balanced OT in Feydy et al. [2019b] and UW in Séjourné et al. [2019]. The design of efficient numerical solvers for CGW is also an interesting avenue for future works, as well as the study of its induced topology.

heading: Acknowledgements
text: The works of Thibault Séjourné and Gabriel Peyré is supported by the ERC grant NORIA. The work of G. Peyré was supported in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d’avenir" program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute).
The authors thank Rémi Flamary for his remarks and advices, as well as Laetitia Chapel for her help to reproduce her experiments.

Review summaries:
The paper extends the Gromov-Wasserstein (GW) distance framework to two novel relaxations, the unbalanced GW (UGW) and the conic GW (CGW), which both address the question of comparing arbitrary positive measures on different metric spaces. The paper also proves that CGW is upper bounded by UGW, and introduces an algorithm to compute UGW based on Sinkhorn iterations. Numerical experiments are presented to illustrate the performance of the new algorithm. 

The paper deals with unbalanced Gromov-Wasserstein (UGW)  to compare un-registered metric measure spaces not necessarily endowed with a probability distribution. It proposes the relaxation of the  usual marginal constraints using quadratic  divergences leading to  a new  divergence that can be used to compare metric measure spaces.  By  adding an adapted entropic regularization to  the UGW formulation, an alternate optimization scheme based on Sinkhorn iterations  is designed. The  paper also  introduces another conic  formulation for Unbalanced GW (CGW) that  defines a distance invariant to isometries. The former is shown  to upper-bound CGW which optimization procedure is computationally  heavy. Both CGW and UGW  are illustrated on simulated. Furthermore application of  UGW to PU  learning shows  the effectiveness  of the  proposed UGW divergence. 

This paper proposes two unbalanced Gromov-Wasserstein formulations, namely Unbalanced Gromov-Wasserstein (UGW) divergence and Conic Gromov-Wasserstein (CGW) distance. Both of them are positive and definite. UGW has a scalable and GPU-friendly algorithm, which is applicable to large learning problems. CGW is a distance between mm-spaces up to isometry. It is theoretically and empirically shown that UGW can be used as a surrogate upper-bounding CGW. The authors then apply UGW to positive unlabeled learning with Caltech office dataset and achieve comparable or better performance than Chapel et al., 2020.