Title: SCALABLE UNBALANCED OPTIMAL TRANSPORT USING GENERATIVE ADVERSARIAL NETWORKS

Abstract: Generative adversarial networks (GANs) are an expressive class of neural generative models with tremendous success in modeling high-dimensional continuous measures. In this paper, we present a scalable method for unbalanced optimal transport (OT) based on the generativeadversarial framework. We formulate unbalanced OT as a problem of simultaneously learning a transport map and a scaling factor that push a source measure to a target measure in a cost-optimal manner. We provide theoretical justification for this formulation, showing that it is closely related to an existing static formulation by Liero et al. (2018) . We then propose an algorithm for solving this problem based on stochastic alternating gradient updates, similar in practice to GANs, and perform numerical experiments demonstrating how this methodology can be applied to population modeling.

1. INTRODUCTION:

We consider the problem of unbalanced optimal transport: given two measures, find a cost-optimal way to transform one measure to the other using a combination of mass variation and transport. Such problems arise, for example, when modeling the transformation of a source population into a target population (Figure 1a ). In this setting, one needs to model mass transport to account for the features that are evolving, as well as local mass variations to allow sub-populations to become more or less prominent in the target population (Schiebinger et al., 2017) .

Classical optimal transport (OT) considers the problem of pushing a source to a target distribution in a way that is optimal with respect to some transport cost without allowing for mass variations. Modern approaches are based on the Kantorovich formulation (Kantorovich, 1942) , which seeks the optimal probabilistic coupling between measures and can be solved using linear programming methods for discrete measures. Recently, Cuturi (2013) showed that regularizing the objective using an entropy term allows the dual problem to be solved more efficiently using the Sinkhorn algorithm. Stochastic methods based on the dual objective have been proposed for the continuous setting (Genevay et al., 2016; Seguy et al., 2017; Arjovsky et al., 2017) . Optimal transport has been applied to many areas, such as computer graphics (Ferradans et al., 2014; Solomon et al., 2015) and domain adaptation (Courty et al., 2014; 2017) .

In many applications where a transport cost is not available, transport maps can also be learned using generative models such as generative adversarial networks (GANs) (Goodfellow et al., 2014) , which push a source distribution to a target distribution by training against an adversary. Numerous transport problems in image translation (Mirza & Osindero, 2014; Zhu et al., 2017; Yi et al., 2017) , natural language translation (He et al., 2016) , domain adaptation (Bousmalis et al., 2017) and biological data integration (Amodio & (a) (b)

Figure 1 : (a) Illustration of the problem of modeling the transformation of a source population µ to a target population ν. In this example, one sub-population is growing more rapidly than the others. (b) Schematic of Monge-like formulations of unbalanced optimal transport. The objective is to learn a transport map T (for transporting mass) and scaling factor ξ (for mass variation) to push the source µ to the target ν, using a deterministic transport map (top) (Chizat et al., 2015) or a stochastic transport map (bottom).

Krishnaswamy, 2018) have been tackled using variants of GANs, with strategies such as conditioning or cycle-consistency employed to enforce correspondence between original and transported samples. However, all these methods conserve mass between the source and target and therefore cannot handle mass variation.

Several formulations have been proposed for extending the theory of OT to the setting where the measures can have unbalanced masses (Chizat et al., 2015; 2018; Kondratyev et al., 2016; Liero et al., 2018; Frogner et al., 2015) . In terms of numerical methods, a class of scaling algorithms (Chizat et al., 2016) that generalize the Sinkhorn algorithm for balanced OT have been developed for approximating the solution to optimal entropy-transport problems; this formulation of unbalanced OT by Liero et al. (2018) corresponds to the Kantorovich OT problem in which the hard marginal constraints are relaxed using divergences to allow for mass variation. In practice, these algorithms have been used to approximate unbalanced transport plans between discrete measures for applications such as computer graphics (Chizat et al., 2016) , tumor growth modeling (Chizat & Di Marino, 2017 ) and computational biology (Schiebinger et al., 2017) . However, while optimal entropy-transport allows mass variation, it cannot explicitly model it, and there are currently no methods that can perform unbalanced OT between continuous measures.

Contributions. Inspired by the recent successes of GANs for high-dimensional transport problems, we present a novel framework for unbalanced optimal transport that directly models mass variation in addition to transport. Concretely, our contributions are the following:

• We propose to solve a Monge-like formulation of unbalanced OT, in which the goal is to learn a stochastic transport map and scaling factor to push a source to a target measure in a cost-optimal manner. This generalizes the unbalanced Monge OT problem by Chizat et al. (2015) .

• By relaxing this problem, we obtain an alternative form of the optimal entropy-transport problem by Liero et al. (2018) , which confers desirable theoretical properties.

• We develop scalable methodology for solving the relaxed problem. Our derivation uses a convex conjugate representation of divergences, resulting in an alternating gradient descent method similar to GANs (Goodfellow et al., 2014) .

• We demonstrate in practice how our methodology can be applied towards population modeling using the MNIST and USPS handwritten digits datasets, the CelebA dataset, and a recent single-cell RNA-seq dataset from zebrafish embrogenesis.

In addition to these main contributions, for completeness we also propose a new scalable method (Algorithm 2) in the Appendix for solving the optimal-entropy transport problem by Liero et al. (2018) in the continuous setting. The algorithm extends the work of Seguy et al. (2017) to unbalanced OT and is a scalable alternative to the algorithm of Chizat et al. (2016) for very large or continuous datasets.

2. PRELIMINARIES:

Notation. Let X , Y ⊆ R n be topological spaces and let B denote the Borel σ-algebra. Let M 1 + (X ), M + (X ) denote respectively the space of probability measures and finite non-negative measures over X . For a measurable function T , let T # denote its pushforward operator: if µ is a measure, then T # µ is the pushforward measure of µ under T . Finally, let π X , π Y be functions that project onto X and Y; for a joint measure γ ∈ M + (X × Y), π X # γ and π Y # γ are its marginals with respect to X and Y respectively. Optimal transport (OT) addresses the problem of transporting between measures in a cost-optimal manner. Monge (1781) formulated this problem as a search over deterministic transport maps. Specifically, given

EQUATION

subject to the constraint T # µ = ν. While the optimal T has an intuitive interpretation as an optimal transport map, the Monge problem is non-convex and not always feasible depending on the choices of µ and ν. The Kantorovich OT problem is a convex relaxation of the Monge problem that formulates OT as a search over probabilistic transport plans. Given

EQUATION

Note that the conditional probability distributions γ y|x specify stochastic maps from X to Y and can be considered a "one-to-many" version of the deterministic map from the Monge problem. In terms of numerical methods, the relaxed problem is a linear program that is always feasible and can be solved in O(n 3 ) time for discrete µ, ν. Cuturi (2013) recently showed that introducing entropic regularization results in a simpler dual optimization problem that can be solved efficiently using the Sinkhorn algorithm. Based on the entropyregularized dual problem, Genevay et al. (2016) and Seguy et al. (2017) proposed stochastic algorithms for computing transport plans that can handle continuous measures.

Unbalanced OT. Several formulations that extend classical OT to handle mass variation have been proposed (Chizat et al., 2015; 2018; Kondratyev et al., 2016) . Existing numerical methods are based on a Kantorovichlike formulation known as optimal-entropy transport (Liero et al., 2018) . This formulation is obtained by relaxing the marginal constraints of (2) using divergences as follows: given two positive measures µ ∈ M + (X ) and ν ∈ M + (Y) and a cost function c : X × Y → R + , optimal entropy-transport finds a measure γ ∈ M + (X × Y) that minimizes

EQUATION

where D ψ1 , D ψ2 are ψ-divergences induced by ψ 1 , ψ 2 . The ψ-divergence between non-negative finite measures P, Q over T ⊆ R d induced by a lower semi-continuous, convex entropy function ψ : R → R ∪ {∞} is

EQUATION

where ψ ∞ := lim s→∞ ψ(s) s and dP dQ Q + P ⊥ is the Lebesgue decomposition of P with respect to Q. Note that mass variation is allowed since the marginals of γ are not constrained to be µ and ν. In terms of numerical methods, the state-of-the-art in the discrete setting is a class of iterative scaling algorithms (Chizat et al., 2016) that generalize the Sinkhorn algorithm for computing regularized OT plans (Cuturi, 2013) . There are no practical algorithms for unbalanced OT between continuous measures, especially in high-dimensional spaces.

3. SCALABLE UNBALANCED OT USING GANS:

In this section, we propose the first algorithm for unbalanced OT that directly models mass variation and can be applied towards transport between high-dimensional continuous measures. The starting point of our development is the following Monge-like formulation of unbalanced OT, in which the goal is to learn a stochastic transport map and scaling factor to push a source to a target measure in a cost-optimal manner.

Unbalanced Monge OT. Let c 1 : X × Y → R + be the cost of transport and c 2 : R + → R + the cost of mass variation. Let the probability space (Z, B(Z), λ) be the source of randomness in the transport map T . Given two positive measures µ ∈ M + (X ) and ν ∈ M + (Y), we seek a transport map T : X × Z → Y and a scaling factor ξ :

EQUATION

subject to the constraint T # (ξµ × λ) = ν. Concretely, the first and second terms of ( 5) penalize the cost of mass transport and variation respectively, and the equality constraint ensures that (T, ξ) pushes µ to ν exactly. A special case of ( 5) is the unbalanced Monge OT problem by Chizat et al. (2015) , which employs a deterministic transport map (Figure 1b ). We consider the more general case of stochastic (i.e. one-to-many) maps because it is a more suitable model for many practical problems. For example, in cell biology, it is natural to think of one cell in a source population as potentially giving rise to multiple cells in a target population.

In practice, one can take Z = R n and λ to be the standard Gaussian measure if a stochastic map is desired; otherwise λ can be set to a deterministic distribution. The following are examples of problems that can be modeled using unbalanced Monge OT. Example 3.1 (Figure 2a ). Suppose the objective is to model the transformation from a source measure (Column 1) to the target measure (Column 2), which represent a population of interest at two distinct time points. The transport map T models the transport/movement of points from the source to the target, while the scaling factor ξ models the growth (replication) or shrinkage (death) of these points. Different models of transformation are optimal depending on the relative costs of mass transport and variation (Columns 3-6). Example 3.2 (Figure 2b ). Suppose the objective is to transport points from a source measure (1st panel, color) to a target measure (1st panel, grey) in the presence of class imbalances. A pure transport map would muddle together points from different classes, while an unbalanced transport map with a scaling factor is able to ameliorate the class imbalance (2nd panel). In this case, the scaling factor tells us explicitly how to downweigh or upweigh samples in the source distribution to balance the classes with the target distribution (3rd panel).

Relaxation. From an optimization standpoint, it is challenging to satisfy the constraint T # (ξµ × λ) = ν. We hence consider the following relaxation of ( 5) using a divergence penalty in place of the equality constraint:

EQUATION

using an appropriate choice of ψ that satisfies the requirements of Lemma C.2 in the Appendix1 .This relaxation is the Monge-like version of the optimal-entropy transport problem (3) by Liero et al. (2018) . Specifically, (T, ξ) specifies a joint measure γ ∈ M + (X × Y) given by

γ(C) := X Z 1 C (x, T (x, z))dλ(z)ξ(x)dµ(x), ∀C ∈ B(X ) × B(Y),

and by reformulating ( 6) in terms of γ instead of (T, ξ), one obtains the objective function for optimal-entropy transport. The main difference between the formulations is their search space, since not all joint measures γ ∈ M + (X × Y) can be specified by some choice of (T, ξ). For example, if T is a deterministic transport map, then γ is necessarily restricted to the set of deterministic couplings. Even if T is sufficiently random, it is generally not possible to specify all joint measures γ ∈ M + (X × Y): in the asymmetric Monge formulation ( 6), all the mass transported to Y must come from somewhere within the support of µ, since the scaling factor ξ allows mass to grow but not to materialize outside of its original support. Therefore equivalence can be established in general only when restricting the support of γ to supp(µ) × Y as described in the following lemma, whose proof is given in the Appendix. Lemma 3.3. Let G be the set of joint measures supported on supp(µ) × Y, and define

Wc,ψ1,ψ2 (µ, ν) := inf γ∈G c dγ + D ψ1 (π X # γ|µ) + D ψ2 (π Y # γ|ν). If (Z, B(Z), λ) is atomless and c 2 is an entropy function, then L ψ (µ, ν) = Wc1,c2,ψ (µ, ν).

Based on the relation between ( 3) and ( 6), several theoretical results for (6) follow from the analysis of optimal entropy-transport by Liero et al. (2018) . Importantly, one can show the following theorem, namely that for an appropriate and sufficiently large choice of divergence penalty, solutions of the relaxed problem (6) converge to solutions of the original problem (5) . The proof is given in the Appendix. Theorem 3.4. Suppose c 1 , c 2 , ψ satisfy the existence assumptions of Proposition B.1 in the Appendix, and let (Z, B(Z), λ) be an atomless probability space. Furthermore, let ψ be uniquely minimized at ψ(1) = 0.

Then for a sequence 0

< ζ 1 < • • • < ζ k < • • • diverging to ∞ indexed by k, lim k→∞ L ζ k ψ (µ, ν) = L(µ, ν).

Additionally, let γ k be the joint measure specified by a minimizer of L ζ k ψ (µ, ν). If L(µ, ν) < ∞, then up to extraction of a subsequence, γ k converges weakly to γ, the joint measure specified by a minimizer of L(µ, ν).

Algorithm. Using the relaxation of unbalanced Monge OT in (6), we now show that the transport map and scaling factor can be learned by stochastic gradient methods. While the divergence term cannot easily be minimized using the definition in (4), we can write it as a penalty witnessed by an adversary function f : Y → (-∞, ψ ∞ ] using the convex conjugate representation (see Lemma B.2):

EQUATION

where ψ * is the convex conjugate of ψ. The objective in (6) can now be optimized using alternating stochastic gradient updates after parameterizing T , ξ, and f with neural networks; see Algorithm 1 2 . The optimization procedure is similar to GAN training and can be interpreted as an adversarial game between (T, ξ) and f :

• T takes a point x ∼ µ and transports it from X to Y by generating T (x, z) where z ∼ λ.

• ξ determines the importance weight of each transported point.

• Their shared objective is to minimize the divergence between transported samples and real samples from ν that is measured by the adversary f . • Additionally, cost functions c 1 and c 2 encourage T, ξ to find the most cost-efficient strategy.

Algorithm 1 Generative-Adversarial Framework for Unbalanced Monge OT Input: Initial parameters θ, φ, ω; step size η; normalized measures μ, ν, constants c µ , c ν .

Output: Updated parameters θ, φ, ω. while (θ, φ, ω) not converged do Sample x 1 , • • • , x n from μ, y 1 , • • • , y n from ν, z 1 , • • • , z n from λ; (θ, φ, ω) := 1 n n i=1 [c µ c 1 (x i , T θ (x i , z i ))ξ φ (x i ) + c µ c 2 (ξ φ (x i )) (8) + c µ ξ φ (x i )f ω (T θ (x i , z i )) -c ν ψ * (f ω (y i )).]

Update ω by gradient descent on -(θ, φ, ω). Update θ, φ by gradient descent (θ, φ, ω). end while Table 1 in the Appendix provides some examples of divergences with corresponding entropy functions and convex conjugates that can be plugged into (7). Further practical considerations for implementation and training are discussed in Appendix C.

Relation to other approaches. The probabilistic Monge-like formulation ( 6) is similar to the Kantorovichlike entropy-transport problem (3) in theory, but they result in quite different numerical methods in practice. Algorithm 1 solves the non-convex formulation (6) and learns a transport map T and scaling factor ξ parameterized by neural networks, enabling scalable optimization using stochastic gradient descent. The networks are immediately useful for many practical applications; for instance, it only requires a single forward pass to compute the transport and scaling of a point from the source domain to the target. Furthermore, the neural architectures of T, ξ imbue their function classes with a particular structure, and when chosen appropriately, enable effective learning of these functions in high-dimensional settings. Due to the non-convexity of the optimization problem, however, Algorithm 1 is not guaranteed to find the global optimum. In contrast, the scaling algorithm of Chizat et al. (2016) based on (3) solves a convex optimization problem and is proven to converge, but is currently only practical for discrete problems and has limited scalability. For completeness, in Section A of the Appendix, we propose a new stochastic method based on the same dual objective as Chizat et al. (2016) that can handle transport between continuous measures (Algorithm 2 in the Appendix). This method generalizes the approach of Seguy et al. (2017) for handling transport between continuous measures and overcomes the scalability limitations of Chizat et al. (2016) . However, the output is in the form of the dual solution, which is less interpretable for practical applications compared to the output of Algorithm 1. In particular, while one can compute a deterministic transport map known as a barycentric projection from the dual solution, it is unclear how best to obtain a scaling factor or a stochastic transport map that can generate samples outside of the target dataset. In the numerical experiments of Section 4, we show the advantage of directly learning a transport map and scaling factor using Algorithm 1.

The problem of learning a scaling factor (or weighting factor) that "balances" measures µ and ν also arises in causal inference. Generally, µ is the distribution of covariates from a control population and ν is the distribution from a treated population. The goal is to scale the importance of different members from the

(a) (b) (c) (d)

Figure 3 : Learning weights on MNIST data using unbalanced OT.

control population based on how likely they are to be present in the treated population, in order to eliminate selection biases in the inference of treatment effects. Kallus (2018) proposed a generative-adversarial method for learning the scaling factor, but they do not consider transport.

4. NUMERICAL EXPERIMENTS:

In this section, we illustrate in practice how Algorithm 1 performs unbalanced OT, with applications geared towards population modeling.

MNIST-to-MNIST. We first apply Algorithm 1 to perform unbalanced optimal transport between two modified MNIST datasets. The source dataset consists of regular MNIST digits with the class distribution shown in column 1 of Figure 3a . The target dataset consists of either regular (for the experiment in Figure 3b ) or dimmed (for the experiment in Figure 3c ) MNIST digits with the class distribution shown in column 2 of Figure 3a . The class imbalance between the source and target datasets imitates a scenerio in which certain classes (digits 0-3) become more popular and others (6-9) become less popular in the target population, while the change in brightness is meant to reflect population drift. We evaluated Algorithm 1 on the problem of transporting the source distribution to the target distribution, enforcing a high cost of transport (w.r.t. Euclidean distance). In both cases, we found that the scaling factor over each of the digit classes roughly reflects its ratio of imbalance between the source and target distributions (Figure 3b-c ). These experiments validate that the scaling factor learned by Algorithm 1 reflects the class imbalances and can be used to model growth or decline of different classes in a population. Figure 3d is a schematic illustrating the reweighting that occurs during unbalanced OT.

MNIST-to-USPS. Next, we apply unbalanced OT from the MNIST dataset to the USPS dataset. As before, these two datasets are meant to imitate a population sampled at two different time points, this time with a large degree of evolution. We use Algorithm 1 to model the evolution of the MNIST distribution to the USPS distribution, taking as transport cost the Euclidean distance between the original and transported images. A summary of the unbalanced transport is visualized in Figure 4a . Each arrow originates from a real MNIST image and points towards the predicted appearance of this image in the USPS dataset. The size of the image reflects the scaling factor of the original MNIST image, i.e. whether it is relatively increasing or decreasing in prominence in the USPS dataset compared to the MNIST dataset according to the unbalanced OT model. Even though the Euclidean distance is not an ideal measure of correspondence between MNIST and USPS digits, many MNIST digits were able to preserve their likeness during the transport (Figure 4b ). We analyzed which MNIST digits were considered as increasing or decreasing in prominence by the model. The MNIST digits with higher scaling factors were generally brighter (Figure 4c ) and covered a larger area of pixels (Figure 4d ) compared to the MNIST digits with lower scaling factors. These results are consistent with the observation that the target USPS digits are generally brighter and contain more pixels. CelebA-Young-to-CelebA-Aged. We applied Algorithm 1 on the CelebA dataset to perform unbalanced OT from the population of young faces to the population of aged faces. This synthetic problem imitates a real application of interest, which is modeling the transformation of a population based on samples taken from two timepoints. Since the Euclidean distance between two faces is a poor measure of semantic similarity, we first train a variational autoencoder (VAE) (Kingma & Welling, 2013) on the CelebA dataset and encode all samples into the latent space. We then apply Algorithm 1 to perform unbalanced OT from the encoded young to the encoded aged faces, taking the transport cost to be the Euclidean distance in the latent space. A summary of the unbalanced transport is visualized in Figure 5a . Each arrow originates from a real face from the young population and points towards the predicted appearance of this face in the aged population. Generally, the transported faces retain the most salient features of the original faces (Figure 5b ), although there are exceptions (e.g. gender swaps) which reflects that some features are not prominent components of the VAE encodings.

Interestingly, the young faces with higher scaling factors were significantly enriched for males compared to young faces with lower scaling factors; 9.6% (9,913/103,287) of young female faces had a high scaling factor as compared to 18.5% (8,029/53,447) for young male faces (Figure 5c , top, p = 0). In other words, our model predicts growth in the prominence of male faces compared to female faces as the CelebA population evolves from young to aged. After observing this phenomenon, we confirmed based on checking the ground truth labels that there was indeed a strong gender imbalance between the young and aged populations: while the young population is predominantly female, the aged population is predominantly male (Figure 5c , bottom).

Zebrafish embroygenesis. A problem of great interest in biology is lineage tracing of cells between different developmental stages or during disease progression. This is a natural application of transport in which the source and target distributions are unbalanced: some cells in the earlier stage are more poised to develop into cells seen in the later stage. To showcase the relevance of learning the scaling factor, we apply Algorithm 1 to recent single-cell gene expression data from two stages of zebrafish embryogenesis (Farrell et al., 2018) . The source population is from a late stage of blastulation and the target population from an early stage of gastrulation (Figure 6a ). The results of the transport are plotted in Figure 6b -c after dimensionality reduction by PCA and T-SNE (Maaten & Hinton, 2008) . To assess the scaling factor, we extracted the cells from the blastula stage with higher scaling factors (i.e. over 90th percentile) and compared them to the remainder of the cells using differential gene expression analysis, producing a ranked list of upregulated genes. Using the GOrilla tool (Eden et al., 2009) , we found that the cells with higher scaling factors were significantly enriched for genes associated with differentiation and development of the mesoderm (Figure 6d ). This experiment shows that analysis of the scaling factor can be applied towards interesting and meaningful biological discovery. Published as a conference paper at ICLR 2019 Algorithm 2 SGD for Unbalanced OT Input: Initial parameters θ, φ; step size η; regularization parameter ; constants c µ , c ν and normalized measures μ, ν Output: Updated parameters θ, φ while (θ, φ) not converged do Sample (x 1 , y 1 ),

• • • , (x n , y n ) from μ ⊗ ν (θ, φ) := 1 n n i=1 [c µ ψ * 1 (-u(x i )) + c ν ψ * 2 (-v(y i )) + c µ c ν e (u(xi)+v(yi)-c(xi,yi))/ ]

Update θ, φ by gradient descent on (θ, φ) end while

The dual solution (u * , v * ) learned from Algorithm 2 can be used to reconstruct the primal solution γ * based on the relation in (10). Concretely, γ * is a transport map that indicates the amount of mass transported between every pair of points in X and Y. Note that the marginals of γ * with respect to X and Y are not necessarily µ and ν, which is where mass variation is implicitly built into the problem. Given γ * , it is possible to also learn an "averaged" deterministic mapping from X to Y. A standard approach is to take the barycentric projection T : X → Y, defined as,

T (x) = min z∈Y E y∼γ * (•|x) d(z, y),

with respect to some distance d : Y × Y → R + . Seguy et al. (2017) proposed a stochastic algorithm for learning such a map from the dual solution, which we reproduce in Algorithm 3.

Algorithm 3 Learning Barycentric Projection:

Input: Learned functions u, v; initial T θ ; distance function d Output:

Updated T θ while T θ not converged do Sample (x 1 , y 1 ), • • • , (x n , y n ) from μ ⊗ ν (θ) := 1 n n i=1 d(x i , T θ (y i ))e (u(xi)+v(yi)-c(xi,yi))/

Update T θ by gradient descent on (θ) end while

B SUPPLEMENT TO SECTION 3 B.1 RELAXATION TO OPTIMAL-ENTROPY TRANSPORT:

The objectives in ( 6) and ( 3) are equivalent if one reformulates (6) in terms of γ instead of (T, ξ), where γ

∈ M + (X × Y) is a joint measure given by γ(C) := X Z 1 C (x, T (x, z))dλ(z)ξ(x)dµ(x), ∀C ∈ B(X ) × B(Y). (12)

Furthermore, the formulations are equivalent if one restricts the search space of (3) to contain only those joint measures that can be specified by some (T, ξ) . This relation between the formulations is formalized by Lemma 3.3.

Proof of Lemma 3.3. First we show L ψ (µ, ν) ≥ Wc1,c2,ψ (µ, ν). If L ψ (µ, ν) = ∞, this is trivial, so assume L ψ (µ, ν) < ∞. Let (T, ξ) be any solution and define γ by ( 12). Note by this definition that

π X # γ(A) = X 1 A (x)ξ(x)dµ(x), ∀A ∈ B(X ), i.e. ξ is the Radon-Nikodym derivative of π X # γ with respect to µ. Also π Y # γ(B) = X Z 1 B (T (x, z))dλ(z)ξ(x)dµ(x) = T # (ξµ × λ)(B), ∀B ∈ B(Y).

It follows that

X Z c 1 (x, T (x, z))dλ(z) ξ(x)dµ(x) + X c 2 (ξ(x))dµ(x) + D ψ (T # (ξµ × λ)|ν) = X ×Y c 1 (x, y)dγ(x, y) + X c 2 (ξ(x))dµ(x) + D ψ (T # (ξµ × λ)|ν)

(by definition of γ, linearity and monotone convergence)

= X ×Y c 1 (x, y)dγ(x, y) + X c 2 ( dπ X # γ dµ )dµ(x) + D ψ (T # (ξµ × λ)|ν) (since ξ is the Radon-Nikodym derivative) = X ×Y c 1 (x, y)dγ(x, y) + D c2 (π X # γ|µ) + D ψ (T # (ξµ × λ)|ν) = X ×Y c 1 (x, y)dγ(x, y) + D c2 (π X # γ|µ) + D ψ (π Y # γ|ν)

≥ Wc,ψ1,ψ2 (µ, ν). Since this inequality holds for any (T, ξ), taking the infimum over the left-hand side yields L ψ (µ, ν) ≥ Wc,ψ1,ψ2 (µ, ν).

To show L ψ (µ, ν) ≤ Wc,ψ1,ψ2 (µ, ν), assume Wc,ψ1,ψ2 (µ, ν) < ∞ and let γ be any solution. By the disintegration theorem, there exists a family of probability measures

{γ y|x } x∈X in M 1 + (Y) such that γ(C) = X Y 1 C (x, y)dγ y|x (y)dπ X # γ, ∀C ∈ B(X ) × B(Y).

Since (Z, B(Z), λ) is atomless, it follows from Proposition 9.1.2 and Theorem 13.1.1 in (Dudley, 2018) that there exists a family of measurable functions {T x : Z → Y} x∈X such that γ y|x is the pushforward measure of λ under T x for all x ∈ X . Denoting T (x, z) : (x, z) → T x (z), then by a change of variables,

γ(C) = X Z 1 C (x, T (x, z))dλ(z)dπ X # γ, ∀C ∈ B(X ) × B(Y).

By hypothesis, π X # γ is restricted to the support of µ, i.e. π X # γ ≪ µ. Let ξ be the Radon-Nikodym derivative dπ X # γ dµ . It follows from the Radon-Nikodym theorem that (T, ξ) satisfy

γ(C) = X Z 1 C (x, T (x, z))dλ(z)ξ(x)dµ(x), ∀C ∈ B(X ) × B(Y),

which is the same relation as in (12) . Same as before,

π X # γ(A) = X 1 A (x)ξ(x)dµ(x), ∀A ∈ B(X ),

and

π Y # γ(B) = X Z 1 B (T (x, z))dλ(z)ξ(x)dµ(x) = T # (ξµ × λ)(B), ∀B ∈ B(Y).

It then follows that

X ×Y c 1 dγ + D c2 (π X # γ|µ) + D ψ (π Y # γ|ν) = X Y c 1 (x, y)dγ y|x (y) dπ X # γ + D c2 (π X # γ|µ) + D ψ (π Y # γ|ν)

(by Fubini's Theorem for Markov kernels)

= X Z c(x, T (x, z))dλ(z) dπ X # γ + D c2 (π X # γ|µ) + D ψ (π Y # γ|ν) (by change of variables) = X Z c(x, T (x, z))dλ(z) dπ X # γ + c 2 ( dπ X # γ dµ )dµ(x) + D ψ (π Y # γ|ν) = X Z c(x, T (x, z))dλ(z) ξ(x)dµ(x) + c 2 (ξ(x))dµ(x) + D ψ (π Y # γ|ν) (by definition of ξ) = X Z c(x, T (x, z))dλ(z) dµ(x) + c 2 (ξ(x))dµ(x) + D ψ (T # (ξµ × λ)|ν) ≥ L ψ (µ, ν).

Since this inequality holds for any γ, this implies that Wc1,c2,ψ (µ, ν) ≥ L ψ (µ, ν), which completes the proof.

For certain cost functions and divergences, it can be shown that L ψ defines a proper metric between positive measures µ and ν, i.e. taking c 2 , ψ to be entropy functions corresponding to the KL-divergence and c 1 = log cos 2 + (d(x, y)), then L ψ (µ, ν) corresponds to the Hellinger-Kantorovich (Liero et al., 2018) or the Wasserstein-Fisher-Rao (Chizat et al., 2018) metric between positive measures µ and ν.

Based on Lemma 3.3, the theoretical analysis of Liero et al. (2018) , and standard results on constrained optimization, it can be shown that for an appropriate and sufficiently large choice of divergence penalty, solutions of the relaxed problem (6) converge to solutions of the original problem (Theorem 3.4).

Proof of Theorem 3.4. Since ζ k ψ(s) converges pointwise to the equality constraint ι = (s), which is 0 for s = 1 and ∞ otherwise, by Lemma 3.9 in (Liero et al., 2018) , we have that lim inf k→∞ Wc1,c2,ζ k ψ (µ, ν) ≥ Wc1,c2,ι= (µ, ν). Additionally, Wc1,c2,ζ k ψ (µ, ν) ≤ Wc1,c2,ι= (µ, ν) for any value of k since for any mini-

mizer γ of Wc1,c2,ι= (µ, ν), it holds that π Y # γ = ν. Hence Wc1,c2,ι= (µ, ν) = c 1 dγ + D c2 (π X # γ|µ) + D ζ k ψ2 (π Y # γ|ν) ≥ Wc1,c2,ζ k ψ (µ, ν),

for all k. Therefore, lim k→∞ Wc1,c2,ζ k ψ (µ, ν) = Wc1,c2,ι= (µ, ν), which then by Lemma 3.3 implies the first part of the proposition.

For the second part, by the hypothesis we have that Wc1,c2,ι= (µ, ν) = C < ∞ and as a consequence Wc1,c2,ζ k ψ (µ, ν) ≤ C for all k. Hence, by Proposition 2.10 in (Liero et al., 2018) , the sequence of minimizers γ k is bounded. If the assumptions of Proposition B.1 are satisfied, then the sequence γ k is equally tight. For assumption (ii) this follows by Proposition 2.10 in (Liero et al., 2018) and for assumption (i) this follows by the Markov inequality: for any λ > 0,

γ k ({(x, y) ∈ X × Y|c 1 (x, y) > λ}) ≤ 1 λ c 1 dγ k ≤ C λ .

Since γ k are bounded and equally tight, by an extension of Prokhorov's theorem (Theorem 2.2 of (Liero et al., 2018) ), there exists a subsequence of γ k that is weakly convergent to some γ. Then by lower semicontinuity, we obtain that

c 1 dγ + D c2 (π X # γ|µ) + lim sup k→∞ D ζ k ψ2 (π Y # γ k |ν) ≤ Wc1,c2,ι= (µ, ν) = C Since D ζ k ψ (π Y # γ k |ν) = ζ k D ψ (π Y # γ k |ν) ≥ 0 and ζ k → ∞, for the left side to be finite, D ψ (π Y # γ k |ν) must converge to 0, so D ψ (π Y

# γ|ν) = 0 by lower semicontinuity. Therefore, γ is a minimizer of Wc1,c2,ι= (µ, ν). By construction of the proof of Lemma 3.3, γ k is equivalent to the product measure induced by minimizers of L ζ k ψ (µ, ν), which implies the second part of the proposition.

B.2 CONVEX CONJUGATE FORM OF DIVERGENCES:

In this section, we present the convex conjugate form of ψ-divergence used to rewrite the main objective as a min-max problem. Lemma B.2. For non-negative finite measures P, Q over T ⊂ R d , it holds that

EQUATION

)

where F is a subset of measurable functions {f : T → (-∞, ψ ∞ ]}. Equality holds if and only if ∃f ∈ F such that (i) the restriction of f to the support of Q belongs to the subdifferential of ψ( dP dQ ), i.e. the Radon-Nikodym derivative of P with respect to Q and (ii) f = ψ ∞ over the support of P ⊥ .

We provide a simple proof of this result. A similar result under stronger assumptions was shown in Nguyen et al. (2008) and used by Nowozin et al. (2016) for generative modeling. A rigorous proof can be found in Liero et al. (2018) .

Proof of Lemma B.2. Note that

D ψ (P |Q) = ψ ∞ P ⊥ (T ) + T ψ dP dQ dQ = ψ ∞ P ⊥ (T ) + T sup ξ∈R {ξ dP dQ -ψ * (ξ)}dQ (by defintion of convex conjugate) = ψ ∞ P ⊥ (T ) + T sup ξ∈(-∞,ψ ∞ ] {ξ dP dQ -ψ * (ξ)}dQ (by Lemma B.3) below = T sup ξ∈(-∞,ψ ∞ ] {ξdP ⊥ + ξ dP dQ dQ -ψ * (ξ)dQ} = sup f ∈F T f dP -ψ * (f )dQ.

By first-order optimality conditions, the optimal f over the support of Q is obtained when dP dQ belongs to the subdifferential of ψ * (f ), or equivalently when f belongs to the subdifferential of ψ( dP dQ ). It is straightforward to see that the optimal f over the support of P ⊥ is equal to ψ ∞ , which completes the proof.

Lemma B.3. If ξ > ψ ∞ := lim s→∞ ψ(s) s , then ψ * (ξ) = ∞. Proof. ψ * (ξ) = sup s∈R sξ -ψ(s) ≥ lim s→∞ s(ξ - ψ(s) s ) = ∞ if ξ > ψ ∞ C PRACTICAL CONSIDERATIONS FOR NUMERICAL EXPERIMENTS

Choice of cost functions. Proposition B.1 gives sufficient conditions on c 1 , c 2 for the problem to be wellposed. In practice, it is often convenient the cost of transport, c 1 , to be some measurement of correspondence between X and Y. For example, we can take c 1 (x, y) to be the Euclidean distance between x and y after mapping them to some common feature space. For the cost of mass adjustment, c 2 , it is generally sensible to choose some convex function that vanishes at 1 (i.e. no mass adjustment) and such that lim x→0 c 2 (x) = lim x→∞ c 2 (x) = ∞ to prevent ξ from becoming too small or too large. Any of the entropy functions shown in Table 1 Table 1 : Table of some common ψ-divergences, associated entropy functions ψ, and convex conjugates ψ * for Algorithm 1, partly adapted from (Nowozin et al., 2016) .

Choice of ψ. In Nowozin et al. (2016) , it was shown that any ψ-divergence could be used to train generative models, i.e. to match a generated distribution P to a true data distribution Q. This is due to Jensen's inequality: for any convex lower semi-continous entropy function ψ, D ψ (P |Q) is uniquely minimized when P = Q, where P, Q are probability measures. However, this does not generally hold when P, Q are not probability measures, as illustrated by the following example.

Example C.1. In the original GAN paper, the discrminative objective, sup f log f (x)dP (x) -log(1 -f (x))dQ(x), corresponds to D ψ (P |Q) with ψ(s) = s log s -(s + 1) log(s + 1) (Nowozin et al., 2016) . If P, Q are probability measures, this divergence is equivalent to the Jensen-Shannon divergence and is minimized when P = Q. If P, Q are non-negative measures with unconstrained total mass, the divergence is minimized when P = ∞ and Q = 0.

When P, Q are not probability measures, we require an additional constraint on ψ to ensure that divergence minimization matchces P to Q: Lemma C.2. Suppose P, Q are non-negative finite measures over T ⊆ R n . If ψ(s) attains a unique minimum at s = 1 with ψ(1) = 0 and ψ ∞ > 0, then D ψ (P |Q) = 0 ⇒ P = Q. Otherwise, then P = Q in general when D ψ (P |Q) is minimized.

Proof. Suppose ψ(s) attains a unique minimum at s = 1 with ψ(1) = 0, ψ ∞ > 0, and P = Q over a region with positive measure. It is straightforward to see by the definition in (4) that D ψ (P |Q) > 0, since at least one of the two terms will be strictly positive. Therefore, the first statement holds. For the second statement, suppose either ψ(s) does not attain a unique minimum at s = 1 or ψ ∞ ≤ 0. If ψ(s) attains a minimum at some s = 1, then taking P = s Q results in a divergence that is equal to or less than P = Q. If ψ ∞ ≤ 0, then letting P = Q + P ⊥ where P ⊥ is a positive measure orthogonal to Q results in a divergence that is equal to or less than P = Q.

Table 1 provides some examples of ψ corresponding to common divergences that can be used for unbalanced OT.

Choice of f . According to Lemma B.2, f should belong to a class of functions that maps from Y to (-∞, ψ ∞ ].

In practice, this can be enforced by parameterizing f using a neural network with a final layer that maps to the correct range, also known as an output activation layer (Nowozin et al., 2016) . Table 1 provides some examples of activation layers that can be used.

Choice of neural architectures. For our experiments in Section 4, we used fully-connected feedforward networks with 3 hidden layers and ReLU activations. For T , the output activation layer was a sigmoid function to map the final pixel brightness to the range (0, 1). For ξ, the output activation layer was a softplus function to map the scaling factor weight to the range (0, ∞).

:

Note that D ψ (T # (ξµ × λ)|ν) = 0 T # (ξµ × λ) = ν ingeneral since the total mass of the transported measure is not constrained; for the relaxation to be valid, ψ(s) should attain a unique minimum at s = 1 (see Lemma C.2).

We assume that one has access to samples from µ, ν, and in the setting where µ, ν are not normalized, then samples to the normalized measures μ, ν as well as the normalization constants cµ, cν . These are reasonable assumptions for practical applications: for example, in a biological assay, one might collect cµ cells from time point 1 and cν cells from time point 2. In this case, the samples are the measurements of each cell and the normalization constants are cµ, cν .

Figures and Tables:

Figure 2: Motivating examples for Unbalanced Monge OT.

Figure 4: Unbalanced Optimal Transport from MNIST to USPS digits.

Figure 5: Unbalanced Optimal Transport from Young to Aged CelebA Faces.

Figure 6: Unbalanced OT on Zebrafish Single-Cell Gene Expression Data. (a) Illustration of the blastula and gastrula stages of zebrafish embryogenesis. (b) T-SNE plot (Maaten & Hinton, 2008) of the unbalanced OT results for a subset of datapoints. The color of the transported points indicates the relative magnitude of the scaling factor (black = high, white = low). (c) Same plot as (b), where we have colored the source points instead of the transported points. (d) GOrilla output of significantly enriched processes (Eden et al., 2009) based on ranked list of enriched genes in cells with high scaling factors (black points from (c)) from a differential gene expression analysis. Processes in the graph are organized from more general (upstream) to more specific (downstream).

are reasonable choices.

Table content: <table><tr><td>Name Kullback-Leibler (KL)</td><td>ψ(s) s log s -s + 1</td><td>D ψ (P |Q) dQ dP -dP + dQ log dP</td><td>ψ  *  (s) e s</td><td>ψ ∞ ∞</td><td>Activation Layer N/A</td></tr><tr><td>Pearson χ 2 Hellinger</td><td>(s -1) 2 √ ( s -1) 2</td><td>( dP dQ -1) 2 dQ ( dP dQ -1) 2 dQ</td><td>s 2 4 + s s 1-s</td><td>∞ 1</td><td>N/A 1 -e s</td></tr><tr><td>Jensen-Shannon</td><td>s log s -(s + 1) log s+1 2</td><td>1 2 D KL (P | P +Q 2 ) + 1 2 D KL (Q| P +Q 2 )</td><td>-log(2 -e s )</td><td>log 2</td><td>log(2) -log(1 + e -s ))</td></tr></table>

Peer review (iclr2019):
 REVIEW The authors propose a novel approach to estimate unbalanced optimal transport between sampled measures that scales well in the dimension and in the number of samples. This formulation is based on a formulation of the entropy-transport problems of Liero et al. where the transport map, growth maps and Lagrangian multipliers are parameterized by neural networks. The effectiveness of the approach is shown on some tasks. This is overall an ingenious contribution that opens a venue for interesting uses of optimal transport tools in learning problems (I can think for instance of transfer learning). As such I think the idea would deserve publication. However, I have some concerns with the way the theory is presented and with the lack of discussions on the theoretical limitations. Also, the theory seems a bit disconnected from the practical set up, and this should be emphasized. These concerns are detailed below. REMARKS ON SECTION 3 I think the theoretical part does not exhibit clearly the relationships with previous literature. The formulation proposed in the paper (6) is not new and consists in solving the optimal entropy-transport problem (2) on the set of product measures gamma that are deterministic, i.e. of the form gamma(x,y) = (id x T)_# (xi mu) for some T:X -> Y and xi : X -> R_+ (here (id x T)(x) =(x,T(x)) ) It is classical in optimal transport to switch between convex/transport plan formulation (easier to study) to non-convex/transport map formulations (easier to interpret). (As a technical note, the support restriction in Lemma 3.2 is automatically satisfied for all feasible plans, for super-linear costs c_2=phi_1). More precisely, since the authors introduce a reference measure lambda on a space Z (these objects are not motivated anywhere, but I guess are used to allow for multivalued transport maps?), they look for plans of the form gamma(x,y) = (pi_x x T)_# (xi mu otime lambda) where (pi_x x T)(x,z) = (x,T(x,z) and "otime" yields product measures) (it is likely that similar connections could be made with the "static" formulations in Chizat et al.). Introduced this way, the relationship to previous literature would have been clearer and the theoretical results are simple consequences of the results in Liero et al., who have characterized when optimal solutions of this form exist. Also this contradicts the remark that the authors make that it is better to model "directly mass variation" as their formulation is essentially equivalent. The paragraph "Relation to Unbalanced OT" is, in my opinion, incomplete. The switch to non-convex formulation introduce many differences to convex approaches that are not mentioned: there is no guarantee that a minimizer can be found, there is a bias introduced by the architecture of the neural network, ... Actually, it is this bias that make the formulation useful in high dimension since it is know that optimal transport suffers from the curse of dimensionality (thus it would be useless to try to solve it exactly in high dimension). I suggest to improve this discussion. OTHER REMARKS A small remark: lemma 3.1 is the convex conjugate formula for the phi-divergence in the first argument. I suggest to call it this way to help the reader connect with concepts he or she already knows. Its rigorous proof (with measurability issues properly dealt with) can be found, for instance, in Liero et al. Theorem 2.7. It follows that the central objective (8) is a Lagrangian saddle-point formulation of the problem of Liero et al., where transport plans, scalings and Lagrange multipliers are parameterized by neural networks. I generally think it is best to make the link with previous work as simple as possible. Also, Appendix C lacks details to understand precisely how the experiments where done. It is written : "In practice, [the correct output range] can be enforced by parameterizing f using a neural network with a final layer that maps to the correct range. In practice, we also found that employing a Lipschitz penalty on f stabilizes training." This triggers two remarks: - (i) how precisely is the correct range enforced? This should be stated. - (ii) a Lipschitz penalty on f yields a class of functions which is very unlikely to have the properties of Lemma 3.1 ; in fact, this amounts to replacing the last term in (6) by a sort of "bounded Lipschitz" distance which has very different property from a f-divergence. This makes the theory of section 3 a bit disconnected from the practice of section 4.