GENERATIVE ADVERSARIAL NETWORKS
Year: 2019

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.

heading: 1 INTRODUCTION
text: 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 &
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.

heading: 2 PRELIMINARIES
text: Notation. Let X ,Y ⊆ Rn 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 µ ∈M 1+(X ), ν ∈M 1+(Y) and a cost function c : X × Y → R+, Monge OT seeks a measurable function T : X → Y minimizing
inf T ∫ X c(x, T (x)) dµ(x) (1)
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 µ ∈ M 1+(X ), ν ∈ M 1+(Y) and a cost function c : X × Y → R+, Kantorovich OT seeks a joint measure γ ∈M 1+(X × Y) subject to πX#γ = µ and πY#γ = ν minimizing
W (µ, ν) := inf γ ∫ X×Y c(x, y) dγ(x, y). (2)
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(n3) 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
Wub(µ, ν) := inf γ ∫ X×Y c(x, y) dγ(x, y) +Dψ1(π X #γ|µ) +Dψ2(πY#γ|ν), (3)
whereDψ1 ,Dψ2 areψ-divergences induced byψ1, ψ2. Theψ-divergence between non-negative finite measures P,Q over T ⊆ Rd induced by a lower semi-continuous, convex entropy function ψ : R→ R ∪ {∞} is
Dψ(P |Q) := ψ′∞P⊥(T ) + ∫ T ψ ( dP dQ ) dQ, (4)
where ψ′∞ := lims→∞ ψ(s) s and dP dQQ+ 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.

heading: 3 SCALABLE UNBALANCED OT USING GANS
text: 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 c1 : X × Y → R+ be the cost of transport and c2 : 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 ξ : X → R+ minimizing
L(µ, ν) := inf T,ξ ∫ X ∫ Z c1(x, T (x, z))dλ(z)ξ(x)dµ(x) + ∫ X c2(ξ(x))dµ(x), (5)
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 = Rn 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:
Lψ(µ, ν) := inf T,ξ ∫ X ∫ Z c1(x, T (x, z))dλ(z)ξ(x)dµ(x) + ∫ X c2(ξ(x))dµ(x) +Dψ(T#(ξµ× λ)|ν), (6) 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 1C(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
W̃c,ψ1,ψ2(µ, ν) := inf γ∈G
∫ c dγ +Dψ1(π X #γ|µ) +Dψ2(πY#γ|ν).
If (Z ,B(Z), λ) is atomless and c2 is an entropy function, then Lψ(µ, ν) = W̃c1,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 c1, c2, ψ 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, limk→∞ 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):
Dψ(T#(ξµ× λ)|ν) = sup f ∫ X ∫ Z f(T (x, z))dλ(z)ξ(x)dµ(x)− ∫ Y ψ∗(f(y))dν(y), (7)
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
1Note that Dψ(T#(ξµ× λ)|ν) = 0 ; T#(ξµ× λ) = ν in general 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).
2We 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ν .
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 c1 and c2 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 x1, · · · , xn from µ̃, y1, · · · , yn from ν̃, z1, · · · , zn from λ;
`(θ, φ, ω) := 1
n n∑ i=1 [cµc1(xi, Tθ(xi, zi))ξφ(xi) + cµc2(ξφ(xi)) (8)
+ cµξφ(xi)fω(Tθ(xi, zi))− cνψ∗(fω(yi)).] 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
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.

heading: 4 NUMERICAL EXPERIMENTS
text: 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.

heading: A DUAL STOCHASTIC METHOD
text: In this section, we present a stochastic method for unbalanced OT based on the regularized dual formulation of (Chizat et al., 2015), which can be considered a natural generalization of Seguy et al. (2017). The dual formulation of (3) is given by
sup u,v − ∫ ψ∗1(−u)dµ− ∫ ψ∗2(−v)dν
subject to u ⊕ v ≤ c, where the supremum is taken over functions u : X → [−ψ′1∞,∞] and v : Y → [−ψ′2∞,∞]. This is a constrained optimization problem that is challenging to solve. A standard technique for making the dual problem unconstrained is to add a strongly convex regularization term to the primal objective (Blondel et al., 2017), such as an entropic regularization term (Cuturi, 2013):
Re(γ) = DψKL(γ|µ⊗ ν)
where > 0. Concretely, this term has a “smoothing” effect on the transport plan, in the sense that it encourages plans with high entropy. By the Fenchel-Rockafellar theorem, the dual of the regularized problem is given by,
W ub(µ, ν) := sup u,v − ∫ ψ∗1(−u)dµ− ∫ ψ∗2(−v)dν − ∫ e(u+v−c)/ d(µ⊗ ν), (9)
where the supremum is taken over functions u : X → [−ψ′1∞,∞] and v : Y → [−ψ′2∞,∞], and the relationship between the primal optimizer γ∗ and dual optimizer (u∗, v∗) is given by
dγ∗ = e(u ∗+v∗−c)/ d(µ⊗ ν). (10)
Next, we rewrite (9) in terms of expectations. 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. Based on these assumptions, we have
W ub(µ, ν) = sup u,v −cµEx∼µ̃ψ∗1(−u(x))− cνEy∼ν̃ψ∗2(−v(y))− cµcνEx∼µ̃,y∼ν̃e(u(x)+v(y)−c(x,y))/ . (11)
If ψ∗1 , ψ ∗ 2 are differentiable, we can parameterize u, v with neural networks uθ, vφ and optimize θ, φ using stochastic gradient descent. This is described in Algorithm 2. Note that this algorithm is a generalization of the algorithm of (Seguy et al., 2017) from classical OT to unbalanced OT. Indeed, taking ψ1, ψ2 to be equality constraints, (9) becomes
sup u,v
∫ u dµ+ ∫ v dν − ∫ e(u+v−c)/ d(µ⊗ ν),
which is the dual of the entropy-regularized classical OT problem.
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 (x1, y1), · · · , (xn, yn) from µ̃⊗ ν̃
`(θ, φ) := 1
n n∑ i=1 [cµψ ∗ 1(−u(xi)) + cνψ∗2(−v(yi)) + 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 Ey∼γ∗(·|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 (x1, y1), · · · , (xn, yn) from µ̃⊗ ν̃
`(θ) := 1
n n∑ i=1 d(xi, Tθ(yi))e (u(xi)+v(yi)−c(xi,yi))/
Update Tθ by gradient descent on `(θ) end while

heading: B SUPPLEMENT TO SECTION 3
text: 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 1C(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ψ(µ, ν) ≥ W̃c1,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 1A(x)ξ(x)dµ(x), ∀A ∈ B(X ),
i.e. ξ is the Radon-Nikodym derivative of πX#γ with respect to µ. Also
πY#γ(B) = ∫ X ∫ Z 1B(T (x, z))dλ(z)ξ(x)dµ(x) = T#(ξµ× λ)(B), ∀B ∈ B(Y).
It follows that∫ X (∫ Z c1(x, T (x, z))dλ(z) ) ξ(x)dµ(x) + ∫ X c2(ξ(x))dµ(x) +Dψ(T#(ξµ× λ)|ν)
= ∫ X×Y c1(x, y)dγ(x, y) + ∫ X c2(ξ(x))dµ(x) +Dψ(T#(ξµ× λ)|ν)
(by definition of γ, linearity and monotone convergence)
= ∫ X×Y c1(x, y)dγ(x, y) + ∫ X c2( dπX#γ dµ )dµ(x) +Dψ(T#(ξµ× λ)|ν)
(since ξ is the Radon-Nikodym derivative)
= ∫ X×Y c1(x, y)dγ(x, y) +Dc2(π X #γ|µ) +Dψ(T#(ξµ× λ)|ν)
= ∫ X×Y c1(x, y)dγ(x, y) +Dc2(π X #γ|µ) +Dψ(πY#γ|ν)
≥ W̃c,ψ1,ψ2(µ, ν). Since this inequality holds for any (T, ξ), taking the infimum over the left-hand side yields Lψ(µ, ν) ≥ W̃c,ψ1,ψ2(µ, ν).
To show Lψ(µ, ν) ≤ W̃c,ψ1,ψ2(µ, ν), assume W̃c,ψ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 1C(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 {Tx : Z → Y}x∈X such that γy|x is the pushforward measure of λ under Tx for all x ∈ X . Denoting T (x, z) : (x, z) 7→ Tx(z), then by a change of variables,
γ(C) = ∫ X ∫ Z 1C(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 1C(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 1A(x)ξ(x)dµ(x), ∀A ∈ B(X ),
and
πY#γ(B) = ∫ X ∫ Z 1B(T (x, z))dλ(z)ξ(x)dµ(x) = T#(ξµ× λ)(B), ∀B ∈ B(Y).
It then follows that∫ X×Y c1dγ +Dc2(π X #γ|µ) +Dψ(πY#γ|ν)
= ∫ X (∫ Y c1(x, y)dγy|x(y) ) dπX#γ +Dc2(π X #γ|µ) +Dψ(πY#γ|ν)
(by Fubini’s Theorem for Markov kernels)
= ∫ X (∫ Z c(x, T (x, z))dλ(z) ) dπX#γ +Dc2(π X #γ|µ) +Dψ(πY#γ|ν)
(by change of variables)
= ∫ X (∫ Z c(x, T (x, z))dλ(z) ) dπX#γ + ∫ c2( dπX#γ dµ )dµ(x) +Dψ(π Y #γ|ν)
= ∫ X (∫ Z c(x, T (x, z))dλ(z) ) ξ(x)dµ(x) + ∫ c2(ξ(x))dµ(x) +Dψ(π Y #γ|ν)
(by definition of ξ)
= ∫ X (∫ Z c(x, T (x, z))dλ(z) ) dµ(x) + ∫ c2(ξ(x))dµ(x) +Dψ(T#(ξµ× λ)|ν) ≥ Lψ(µ, ν).
Since this inequality holds for any γ, this implies that W̃c1,c2,ψ(µ, ν) ≥ Lψ(µ, ν), which completes the proof.
Due to the near equivalence of the formulations, several theoretical results for (6) follow from the analysis of optimal entropy-transport by Liero et al. (2018), such as the following existence and uniqueness result:
Proposition B.1. Suppose Lψ(µ, ν) < ∞, c2 is convex and lower semi-continuous, and (Z ,B(Z), λ) is atomless. If (i) c1 has compact sublevel sets in X × Y and c′2∞ + ψ′∞ > 0, or (ii) c′2∞ = ψ′∞ =∞ then a minimizer of Lψ(µ, ν) exists. If c2, ψ are strictly convex, ψ′∞ =∞ and c1 satisfies Corollary 3.6 (Liero et al., 2018), then the joint measure γ specified by any minimizer of Lψ(µ, ν) is unique.
Proof of Proposition B.1. Note that W̃c1,c2,ψ(µ, ν) is equivalent to Wc1,c2,ψ(µ, ν) when X is restricted to the support of µ. If c1, c2, ψ satisfy (i) or (ii), they also satisfy (i) and (ii) when X is restricted to the support of µ. By Theorem 3.3 of (Liero et al., 2018), W̃c1,c2,ψ(µ, ν) has a minimizer. It follows from the construction of the proof of Lemma 3.3 that a minimizer of Lψ(µ, ν) also exists. For uniqueness, if ψ′∞ = ∞, then it follows from Lemma 3.5 of (Liero et al., 2018) and the fact that minimizers are restricted to G that the marginals πX#γ, πY#γ are uniquely determined for any solution γ of W̃c1,c2,ψ(µ, ν). The uniqueness of γ then follows from the proof of Corollary 3.6 in (Liero et al., 2018). It follows from the construction of the proof of Lemma 3.3 that the product measure generated by the minimizers of Lψ(µ, ν) is unique, 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 c2, ψ to be entropy functions corresponding to the KL-divergence and c1 = 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 (5) (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 infk→∞ W̃c1,c2,ζkψ(µ, ν) ≥ W̃c1,c2,ι=(µ, ν). Additionally, W̃c1,c2,ζkψ(µ, ν) ≤ W̃c1,c2,ι=(µ, ν) for any value of k since for any minimizer γ of W̃c1,c2,ι=(µ, ν), it holds that π Y #γ = ν. Hence
W̃c1,c2,ι=(µ, ν) = ∫ c1dγ +Dc2(π X #γ|µ) +Dζkψ2(π Y #γ|ν) ≥ W̃c1,c2,ζkψ(µ, ν),
for all k. Therefore, limk→∞ W̃c1,c2,ζkψ(µ, ν) = W̃c1,c2,ι=(µ, ν), which then by Lemma 3.3 implies the first part of the proposition.
For the second part, by the hypothesis we have that W̃c1,c2,ι=(µ, ν) = C < ∞ and as a consequence W̃c1,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 |c1(x, y) > λ}) ≤ 1
λ
∫ c1dγ
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 ∫
c1dγ̄ +Dc2(π X #γ̄|µ) + lim sup
k→∞ Dζkψ2(π
Y #γ k|ν) ≤ W̃c1,c2,ι=(µ, ν) = C
Since Dζkψ(π Y #γ k|ν) = ζkDψ(π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 W̃c1,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 ⊂ Rd, it holds that
Dψ(P |Q) ≥ sup f∈F
∫ fdP − ψ∗(f)dQ (13)
whereF 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 ψ( dPdQ ), 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 fdP − ψ∗(f)dQ.
By first-order optimality conditions, the optimal f over the support of Q is obtained when dPdQ belongs to the subdifferential of ψ∗(f), or equivalently when f belongs to the subdifferential of ψ( dPdQ ). 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 ξ > ψ′∞ := lims→∞ ψ(s) s , then ψ ∗(ξ) =∞.
Proof.
ψ∗(ξ) = sup s∈R sξ − ψ(s)
≥ lim s→∞ s(ξ − ψ(s) s ) =∞ if ξ > ψ′∞

heading: C PRACTICAL CONSIDERATIONS FOR NUMERICAL EXPERIMENTS
text: Choice of cost functions. Proposition B.1 gives sufficient conditions on c1, c2 for the problem to be wellposed. In practice, it is often convenient the cost of transport, c1, to be some measurement of correspondence between X and Y . For example, we can take c1(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, c2, it is generally sensible to choose some convex function that vanishes at 1 (i.e. no mass adjustment) and such that limx→0 c2(x) = limx→∞ c2(x) =∞ to prevent ξ from becoming too small or too large. Any of the entropy functions shown in Table 1 are reasonable choices.
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 ⊆ Rn. If ψ(s) attains a unique minimum at s = 1 with ψ(1) = 0 and ψ′∞ > 0, then Dψ(P |Q) = 0 ⇒ P = Q. Otherwise, then P 6= 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 6= 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′ 6= 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,∞).
Gradient penalties. The training of GANs using alternating stochastic gradient descent is not guaranteed to converge locally (Mescheder et al., 2018). Stability can be improved using regularization in the form of added instance noise (Arjovsky & Bottou, 2017) or gradient penalties (Gulrajani et al., 2017). We observed that enforcing gradient penalities on ξ and f , while not necessary, improved the stability of Algorithm 1. A gradient penalty on ξ effectively restricts the search to sufficiently smooth ξ, which is reasonable in practice considering that similar regions of the source measure often have similar rates of growth or contraction. A gradient penalty on f changes the nature of the relaxation from (5) to (6): the right-hand side of (7) is no longer equivalent to the ψ-divergence, but is rather a lower-bound with a relation to bounded Lipschitz metrics (Gulrajani et al., 2017). In this case, while the problem formulation is not equivalent to optimal entropy-transport, it is still a valid relaxation of unbalanced Monge OT in (5).
Choice of λ. One can take λ to be the standard Gaussian measure if a stochastic mapping is desired, similar to Almahairi et al. (2018). If a deterministic mapping is desired, then λ is set to a deterministic distribution.
Improved training dynamics. Recall that the objective function for our alternating gradient updates is
`(θ, φ, ω) := 1
n n∑ i=1 [cµc1(xi, Tθ(xi, zi))ξφ(xi) + cµc2(ξφ(xi))
+ cµξφ(xi)fω(Tθ(xi, zi))− cνψ∗(fω(yi))].
Early in training, ξφ can become very small for some xi as none of the transported samples resemble samples from the target distribution. As a result, Tθ may improve very slowly for some inputs xi. One way to address this issue without changing the fixed point is to update ξφ, fω using the above objective and update Tθ using the following objective.
`(θ, φ, ω) := 1
n n∑ i=1 cµ[c1(xi, Tθ(xi, zi)) + fω(Tθ(xi, zi))]
Note that we have omitted terms in the original objective that do not include Tθ, and which therefore do not affect the gradient update. For the terms that remain, the difference is that we are rescaling the contribution of each sample xi by 1/ξ(xi). As long as ξ(xi) > 0, this has the effect of rescaling the gradient update of the loss function with respect to each Tθ(xi, zi) without changing the direction of the update.

Review summaries:
In this paper the authors consider the unbalanced optimal transport problem between two measures with different total mass. The authors introduce first the now standard Kantorovich-like formulation, which considers a coupling whose marginals are penalized to look like the two target measures. The authors introduce a second formulation in (2), somewhat a Kantorovich/Monge hybrid that involves a "random" Monge map where the target point T(x) of a point x now depends also on an additional random variable z, to desribe T(x,z). The authors also consider a local mass creation term (\xi) to weight the initial measure \mu.

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.