TENDED SEMI-DISCRETE OPTIMAL TRANSPORT
Year: 2020

Abstract:
Current generative models like generative adversarial networks (GANs) and variational autoencoders (VAEs) have attracted huge attention due to its capability to generate visual realistic images. However, most of the existing models suffer from the mode collapse or mode mixture problems. In this work, we give a theoretic explanation of the both problems by Figalli’s regularity theory of optimal transportation maps. Basically, the generator compute the transportation maps between the white noise distributions and the data distributions, which are in general discontinuous. However, deep neural networks (DNNs) can only represent continuous maps. This intrinsic conflict induces mode collapse and mode mixture. In order to tackle the both problems, we explicitly separate the manifold embedding and the optimal transportation; the first part is carried out using an autoencoder (AE) to map the images onto the latent space; the second part is accomplished using a GPU-based convex optimization to find the discontinuous transportation maps. Composing the extended optimal transport (OT) map and the decoder, we can finally generate new images from the white noise. This AE-OT model avoids representing discontinuous maps by DNNs, therefore effectively prevents mode collapse and mode mixture.

heading: 1 INTRODUCTION
text: Generative adversarial networks (GANs) (Goodfellow et al. (2014)) and variational autoencoders (VAEs) (Kingma & Welling (2013)) emerge as the dominant approaches for unconditional image generation. When trained on appropriate datasets, they are able to produce realistic and visual appealing samples. GAN methods train an unconditional generator that regresses real images from random noise and a discriminator that measures the difference between generated samples and real images. Despite GANs’ advantages, they have critical drawbacks. 1) Training of GANs are tricky and sensitive to hyperparameters. 2) GANs suffer from mode collapse, in which the generator only learns to generate few modes of data distribution while missing others, although samples from the missing modes occur throughout the training data (see e.g. Goodfellow (2016)). While for the VAEs, the encoder is used to map the data distribution to a Gaussian latent distribution, which is then mapped back to the data distribution by the decoder. While standard VAEs tend to capture all modes, they often generate ambiguous images on multimodal real data distributions. We propose that these phenomena relates deeply with the singularities of distribution transport maps.
Manifold Distribution Hypothesis In deep learning, the manifold distribution hypothesis is well accepted, which assumes the distribution of a specific class of natural data is concentrated on a low dimensional manifold embedded in the high dimensional data space Tenenbaum et al. (2000). Therefore, GANs and VAEs implicitly aim to accomplish two major tasks: 1) manifold embedding: to find the encoding/decoding maps between the data manifold embedded in the image space and the latent space; 2) probability distribution transport: to transport a given white noise distribution to the data distribution, either in the latent or in the image space.
Distribution Transformation The generator of GAN model and the decoder of VAE model are trained to compute a transport map that transforms a known continuous distribution (e.g. Gaussian white noise) to the real data distribution. Namely, the transport map pushes forward the given noise distribution to a generated distribution to approximate the real data distribution, the similarity between the two distributions determines the generalization ability of the generator Ben-David et al. (2010).
Discontinuity and Mode Collapse/Mixture It is a common practice among GAN/VAE models that the generators/decoders are expressed by deep neural networks, which can only represent continuous mappings. Unfortunately, as pointed out by works Nagarajan & Kolter (2017); Khayatkhoei et al. (2018); Xiao et al. (2018), the transport maps may be discontinuous when there are multiple modes in the data distribution. This intrinsic conflict can cause mode collapse or mode mixture. The later means that the generated samples are mixtures of multiple modes and look spurious or ambiguous. 1
Furthermore, even when the real data distribution has a single mode, ambiguous data (e.g. a human face image with mismatched eye colors) can still present. This can be explained by Brenier’s polar factorization theorem Brenier (1991b; 1987; 1991a) and Figalli’s regularity theorem Figalli (2010); Chen & Figalli (2017) (Thm. 5 in Appendix B), which asserts that if the support of the target distribution is not convex, then there will be singularity sets on the support of the source distribution, such that the transport map is discontinuous on these sets. This shows the intrinsic difficulties of conventional GANs/VAEs cannot be eliminated, as shown in Fig. 1.
Conquering Mode Collapse/Mixture However, according to Brenier (1987; 1991a) theorem, the optimal transport map can be represented as the gradient map of the Brenier potential. At the regular points, the Brenier potential is differentiable, its gradient map (the transport map) is continuous; at the singularities, the Brenier potential is continuous but not differentiable, and its gradient map is discontinuous. Conventional GANs and VAEs model the gradient map directly and encounter the trouble of discontinuity. In contrast, we propose to model the globally continuous Brenier potential to avoid mode collapse/mixture.
More specifically, our proposed AE-OT model separates the manifold embedding step and the probability distribution transformation step, the former is carried out by an autoencoder (AE), the latter is accomplished by a convex optimization framework (OT). The OT step computes the Brenier potential explicitly and is able to locate the singularity set (the discontinuous points of the gradient map) based on Figalli’s theory. Our experimental results demonstrate that the proposed method can not only cover all of the modes, but also avoid generating spurious samples (mode mixture).
Contributions (i) From theoretical aspect, this work gives a thorough explanation of mode collapse and mode mixture by the regularity theory of optimal transportation developed by Figalli (2018
1For example, a generator generates obscure digits mixing 0 and 8 but neither 0 nor 8 on the MNIST dataset.
Fields medalist) and the reasons why standard GANs/VAEs cannot solve this problem perfectly. (ii) From practical aspect, this work proposes a novel model called AE-OT, which first encodes the data manifold into the latent space, then compute the Brenier potential to represent the optimal transportation map in the latent space. The Figalli’s singularity set can be located efficiently and avoided when generating new samples. In this way, our model eliminates mode collapse and mode mixture successfully. (iii) The algorithm for finding the Brenier potential and the optimal transportation map can be accelerated with GPU based convex optimization algorithm. The method converges to the unique global optimum with bounded error estimate. (iv) Our experiment results demonstrate the efficiency and efficacy of the proposed method.

heading: 2 RELATED WORK
text: Optimal Transport Optimal transport (OT) has been successfully applied in the areas like economy, optics and machine learning, as surveyed in Solomon (2018) and Peyré & Cuturi (2018). In Gu et al. (2016) , the intrinsic connection between Brenier theory in OT and Alexandroff theory in convex geometry was established, and applied for deep learning in Lei et al. (2017) by an convex optimization. Figalli and the collaborators Figalli (2010); Chen & Figalli (2017) proposed that when the support of the data distribution is non-convex, the transport map will be discontinuous.
Generative models Generative model is one of the main tasks in the machine learning field. One of the most used image generation methods is Variational Autoencoders (VAEs), where the decoders approximate real data distributions from a Gaussian distribution in a variational approach (e.g Kingma & Welling (2013) and Rezende et al. (2014)). Later, Adversarial Autoencoders (AAEs) Makhzani et al. (2015) and Wasserstein Autoencoders (WAEs) Tolstikhin et al. (2018) are proposed following the similar scheme. Although VAEs are relatively simple to train, images they generate look blurry. Generative Adversarial Networks (GANs) Goodfellow et al. (2014) proposed by Goodfellow et.al can produce better quality images. While being a powerful tool in generating realistically looking samples, GANs suffer from the mode collapse problem. To solve this problem, a huge number of methods, including changing the loss function (e.g. Wasserstein GAN Arjovsky et al. (2017)), regularizing the discriminators to be Lipschitz (clipping Arjovsky et al. (2017), gradient regularization Gulrajani et al. (2017), Mescheder et al. (2018) or spectral normalization Miyato et al. (2018)), were proposed.
Besides, various non-adversarial methods has also been proposed recently. GLO Bojanowski et al. (2017) adopts the “encoder-less autoencoder" method to generate new images with a non-adversarial loss function. IMLE Li & Malik (2018) proposed an ICP related generative model training approach. Later GLANN Hoshen & Malik (2019) combines the advantages of GLO and IMLE, an embedding from image space to latent space was first found using GLO and then a transformation between an arbitrary distribution and latent code was computed using IMLE.
Mitigating Mode Collapsing Recently, Nagarajan & Kolter (2017); Khayatkhoei et al. (2018); Xiao et al. (2018) also realize the training difficulties of GANs come from the approximation of discontinuous functions with continuous DNNs. By the gradient-based regularization, GDGAN Nagarajan & Kolter (2017) do relieve the mode collapse phenomenon of GANs, but mode mixture still exists. Khayatkhoei et al. (2018) proposes to use multiple GANs to overcome the mode collapse. Xiao et al. (2018) proposed to embed the images into a latent space according to Bourgain’s theorem, and train the generator by sampling a Gaussian mixture distribution in the latent space instead of a unimodal Gaussian. The recently introduced normalized diversification by Liu et al. (2018) can also help overcome mode collapse successfully. However, all of them cannot solve the mode mixture well.
All these works Nagarajan & Kolter (2017); Khayatkhoei et al. (2018); Xiao et al. (2018) explain that if the target data distribution has multiple modes, the transport map is discontinuous, but DNNs can only represent continuous mappings, the intrinsic conflict causes mode collapse.

heading: 3 COMPUTATIONAL ALGORITHMS
text: Overview of AE-OT Model Our AE-OT model is summarized in Fig. 2, it has two major components: i) (AE) An autoencoder is trained to encode (fθ) the data manifold from the image space X to the latent space Z , and map the data distribution to the latent code distribution; then the decoder gξ decodes the latent code back to the data manifold. ii) (OT) This module computes the optimal transportation map T from the noise distribution to the latent code distribution. First, the Brenier potential is found by a convex optimization process according to Gu et al. (2016), whose gradient is the semi-discrete optimal transport map, where the target is the discrete set of latent codes of training
samples; then the transport map is piece-wise linearly extended to a global continuous map T̃ , where the image domain becomes a simplicial complex obtained by triangulating the above latent codes. Finally, the singularity set in the source domain is located and avoided when generating new samples. As a result, given a random noise x, we can get the generated image by gξ ◦ T̃ (x). Semi-Discrete OT Map Suppose the source measure µ (Gaussian or uniform distribution) is absolutely continuous defined on a convex domain Ω ⊂ Rd, the target domain is a discrete set, Y = {y1, y2, · · · , yn}, yi ∈ Rd, the target measure is a Dirac measure, ν = ∑n i=1 νiδ(y − yi), i =
1, 2, . . . , n, with the equal total mass as the source measure, µ(Ω) = ∑n i=1 νi. Under a semi-
discrete transport map T : Ω → Y , a cell decomposition is induced Ω = ⋃n i=1Wi, such that every x in each cell Wi is mapped to the target yi, T : x ∈ Wi 7→ yi. The map T is measure preserving, denoted as T#µ = ν, if the µ-volume of each cell Wi equals to the ν-measure of the image T (Wi) = yi, µ(Wi) = νi. The cost function is given by c : Ω × Y → R, where c(x, y) represent the cost for transporting a unit mass from x to y. The total cost of T is given by ∫
Ω c(x, T (x))dµ(x) = ∑n i=1 ∫ Wi c(x, yi)dµ(x). Semi-discrete optimal transport map is the
measure-preserving map that minimizes the total cost, T ∗ := arg minT#µ=ν ∫ Ω c(x, T (x))dµ(x).
When the cost function is the L2 distance c(x, y) = 1/2‖x − y‖2, Brenier’s theorem claims that the semi-discrete OT map is given by the gradient map of a piece-wise (PL) convex function, the so-called Brenier potential uh : Ω→ R, uh(x) := maxni=1{πh,i(x)}, where πh,i(x) = 〈x, yi〉+ hi is the hyperplane corresponding to yi ∈ Y . As shown in Fig. 3(a), the projection of the graph of uh decomposes Ω into cells Wi(h), each cell Wi(h) is the projection of the supporting plane πh,i(x). The height vector h is the unique optimizer of the following convex energy under the condition that∑
i hi = 0,
E(h) = ∫ h 0 n∑ i=1 wi(η)dηi − n∑ i=1 hiνi, (1)
where wi(η) is the µ-volume ofWi(η). The convex energyE(h) can be optimized simply by gradient descend method with∇E(h) = (wi(h)− νi)T . The key is to compute the µ-volume wi(h) of each cell Wi(h), which can be estimated using conventional Monte Carlo method. We draw N random samples from µ distribution, {xj} ∼i.i.d. µ, ∀j ∈ J , the estimated µ-volume of each cell is ŵi(h) = #{j ∈ J | xj ∈ Wi(h)}/N . Given xj , we can find Wi in which xj ∈ Wi by i = arg maxi{〈xj , yi〉 + hi}, i = 1, 2, . . . , n. When N is large enough, ŵi(h) converges to wi(h). Then the gradient of the energy is approximated as ∇E ≈ (ŵi(h)− νi)T . Once the gradient is estimated, we can use Adam algorithm Kingma & Ba (2015) to minimize the energy. Sampling of x is independent of each other and finding the cell that
x is located only involves matrix multiplication and sorting. Hence the Monte Carlo method has a natural parallel computation implementation on GPUs.
The approximation error is proportional to the inverse of the square root the amount of Monte Carlo samples. Asymptotically the number of Monte Carlo samples increases exponentially with respect to the dimension d (see e.g. Weed & Bach (2017)). This brings huge computational burdens. To find a good balance between precision and speed, we adaptively adjust the number of random samples. In practice, we apply the following strategy: if the energy E(h) ceases decreasing for a number of consecutive steps, we double the amount of Monte Carlo samples. The algorithmic details of semi-discrete OT map are summarized in Alg. 1.
Algorithm 1 Semi-discrete OT Map 1: Input: Latent codes Y = {yi}i∈I , empirical latent code dis-
tribution ν = 1|I| ∑
i∈I δyi , number of Monte Carlo samples N , positive integer s.
2: Output: Optimal transport map T (·). 3: Initialize h = (h1, h2, . . . , h|I|)← (0, 0, . . . , 0). 4: repeat 5: GenerateN uniformly distributed samples {xj}Nj=1. 6: Calculate∇h = (ŵi(h)− νi)T . 7: ∇h = ∇h−mean(∇h). 8: Update h by Adam algorithm with β1 = 0.9, β2 = 0.5. 9: ifE(h) has not decreased for s steps then
10: N ← N × 2. 11: end if 12: until Converge 13: OT map T (·)← ∇(maxi〈·, yi〉+ hi).
Algorithm 2 Generate latent code 1: Input: Optimal transport map T (·), number of samples to gen-
erate n, angle threshold θ̂. 2: Output: Generated latent code P . 3: Compute ĉi by Monte Carlo method. 4: repeat 5: Sample x ∼ µ, Find the smallest d + 1 vertex around x as
{d(x, ĉi0 ), d(x, ĉi1 ), . . . , d(x, ĉid )}. 6: Compute dihedral angles θik between πi0 and πik . 7: Select θik with θik ≤ θ̂, result in îk = 0, 1, . . . , d1. 8: if ∀k, θik > θ̂ then Abandon x 9: else Generate latent code T̃ (x) = ∑d1 k=0 λkT (ĉîk ) with
λk = d −1(x, ĉîk
)/ ∑d1
j=0 d −1(x, ĉîj ).
10: end if 11: until Generate n new latent code
Piece-wise Linear Extension The semi-discrete OT map ∇uh : Ω → Y maps all x ∈ Ω to the latent codes of training samples {yi}’s and won’t generate new samples. Therefore, we extend the semi-discrete OT map T = ∇uh to a piecewise linear (PL) mapping T̃ as follows. The projection of uh in the source domain induces a cell decomposition of Ω, of which each cell is of µ-volume νi and is mapped to the corresponding yi. By representing the cells by their µ-mass centers as ci := ∫ Wi(h)
xdµ(x), we can get the point-wise map t : ci 7→ yi. The Poincaré of the cell decomposition induces a triangulation of the centers C = {ci}: if Wi ∩Wj 6= ∅, then ci is connected with cj to form an edge [ci, cj ]. Similarly, if Wi0 ∩Wi1 · · · ∩Wik 6= ∅, then there is a k-dimensional simplex [ci0 , ci1 , . . . , cik ]. All these simplices form a triangulation of C (a simplicial complex), denoted as T (C) (the green triangles in the left of Fig. 3(b)). We can triangulate Y in the same way to obtain the triangulation T (Y ) (the green triangles in the right of Fig. 3(b)). Once a random sample x is drawn from the distribution µ, we can find the simplex σ in T (C) containing x. Assume the simplex σ has d+ 1 vertices {ci0 , ci1 , . . . , cid}, the bary-centric coordinates of x in σ is defined as x = ∑d k=0 λkcik , and ∑d k=0 λk = 1 with all λk non-negative. Then the generated latent code of x
under this piece-wise linear map is given by T̃ (x) = ∑d k=0 λkyik (In Fig. 3(b), the green dot x1 is
mapped to be T̃ (x1)). Because all of the yis are used to construct the simplicial complex T (Y ) in the support of the target distribution, we can guarantee that no mode is lost.
In practice, the µ-mass center ci is approximated by the mean value of all the Monte-Carlo samples inside Wi(h), ĉi = ∑ xj∈Wi xj/#{xj ∈ Wi}, where xj ∼ µ. The connectivity information T (C) is too complicated to construct and to store in high dimensional space, thus T (C) is not explicitly built. Instead, we find the simplex σ ∈ T (C) containing x as follows: given a random point x ∈ Ω, evaluate and sort its Euclidean distances to the centers d(x, ĉi), i = 1, 2, . . . , n in the ascending order. Suppose the first d+ 1 items are {d(x, ĉi0), d(x, ĉi1), . . . , d(x, ĉid)}, then σ is formed by {ĉik}. The bary-centric coordinates λ̂ik are estimated as λ̂ik = d −1(x, ĉik)/ ∑d k=0 d
−1(x, ĉik). However, this may generate some spurious samples. To overcome it, we need further to detect the singular set.
Singular Set Detection According to Figalli’s theory Figalli (2010); Chen & Figalli (2017), if there are multiple modes or the support of the target distribution is concave, there will be singular sets Σ ⊂ Ω, where the Brenier potential is continuous but not differentiable, making its gradient map, i.e. the transport map, discontinuous.
As shown in Fig. 3(a), the source distribution is uniformly defined on Ω, and the target empirical distribution has two modes. There is one ridge (the red line) on the Brenier potential uh, whose projection is the singular set Σ (the red line in Ω). Ω \ Σ consists of two connected components, each of them is mapped onto a single mode. Σ consists of codimension 1 facets of cells. If Wi(h) ∩ Wj(h) ⊂ Σ, then the dihedral angle between two supporting planes πh,i and πh,j of uh is prominently large. Therefore, on the graph of Brenier potential, we pick the pairs of facets whose dihedral angles are larger than a given threshold, the projection of their intersection gives a co-dimension 1 cell in the singular set Σ. During the generation process, if a random sample x is around Σ, it will be mapped by T̃ to the gaps among the modes. When generating new latent codes, we just abandon such samples. And this help our method prevent the mode mixture phenomenon.
Given the extended OT map T̃ (x), some of the polyhedrons transverse the singular set (the red lines of Fig .3(b)), which means that different vertices of the polyhedron belongs to different mode. If the sample x falls into such a polyhedron (the dotted red triangle), we just abandon it (as shown in Fig. 3(b), the red dot x2 is just abandoned). Specifically, given x, we can detect if it belongs to the singular set by checking the angles θik between πi0 and πik , k = 1, 2, . . . , d as θik = 〈yi0 , yik〉/‖yi0‖ · ‖yik‖. If all of the angles θik is larger than a threshold θ̂, we say x belongs to the singular set and just abandon it. Or we just select a subset {πik} with θik ≤ θ̂, denoted as {πîk , k = 0, 1, . . . , d1}. Then we can compute λk = d−1(x, ĉîk)/ ∑d1 j=0 d −1(x, ĉîj ) and T̃ (x) = ∑d1 k=0 λkT (ĉîk). Intuitively, T̃ (·) smooths the discrete function T (·) in regions where latent codes are dense and keep the discontinuity of T (·) where latent codes are very sparse. In this way we avoid generating spurious latent code and thus improve the generation quality. The algorithm to generate new code is shown in Alg. 2 and the effect of threshold filtering is further investigated in Appendix 4.1.

heading: 4 EXPERIMENTS
text: In order to validate that the proposed method can solve the mode collapse/mixture problems and generate controllable high quality images, several experiments are conducted.
The first experiment explores the influence of the angle thresholds for the singularity detection on synthetic dataset.
The second experiment focuses on toy sets, so that the complexity of the tasks can be manually controlled and the mode and quality of the generated samples can be accurately computed. Lin et al. (2018) did a large-scale comparison with previous methods that explicitly proposed to mitigate mode collapse and thus established a baseline for comparison. For consistent evaluation, we set up our experiment on the same benchmark dataset as theirs, and make the comparison.
In the last experiment, we run the proposed method mainly on 4 public datasets, MNIST LeCun & Cortes (2010), MNIST-FANSION Han Xiao & Vollgraf (2017), CIFAR-10 Krizhevsky (2009) and CelebA Zhang et al. (2018), just like the authors of Hoshen & Malik (2019) Sajjadi et al. (2018) Lucic et al. (2018) did in their papers. Besides, the architecture of the decoder is the same as Lucic et al. (2018), in which the authors did a large-scale study to evaluated the best performance of 8
different generative models including various GAN models and VAE, and the encoder is set to be the mirror of decoder.

heading: 4.1 SINGLE PARAMETER SELECTIVE INTERPOLATION
text: On synthetic datasets, effects of angle threshold filtering can be visually inspected. As illustrated in Fig. 4, number of mode is a monotonically increasing function with respect to angle threshold θ̂. Quality of generated samples is effected directly by choosing different θ̂. Generally, small θ̂ encourages interpolation in between closely related real samples while too large θ̂ will result in interpolation between samples from different modes, which might in turn lower generation quality. On synthetic datasets, where modes are isotropic and different modes are clearly separable, an ideal θ̂ that captures all modes while avoids generating low quality samples can be chosen within a relatively wide band. For real world datasets of unknown modes, an ideal θ̂ needs to be hand tuned as the separability of different modes depends largely on input data pattern and quality of the embedding map.

heading: 4.2 MITIGATION OF MODE COLLAPSE AND MODE MIXTURE IN SYNTHETIC DATASET
text: Since synthetic dataset consists of explicit distributions and known modes, mode collapse and the quality of the generated sample can be accurately measured. We choose the same synthetic datasets and metrics as in Lin et al. (2018). Specifically, we use 2D-ring and 2D-grid for test sets and Number of modes, Percentage of high-quality samples, reverse Kullback-Leibler (KL) divergence as evaluation metrics. Number of modes counts the amount of modes captured by samples produced a generative model. Percentage of high-quality samples measures the proportion of samples that generated within three standard deviations of the nearest mode. reverse KL divergence measures how well generated samples balance among all modes regarding the real distribution. In Lin et al. (2018), the authors evaluated GAN, ALI, MD and PacGAN on synthetic sets with above three metrics. Each networks are trained under the same generator architecture with a total of approximated 400K training parameters. For AE-OT test, since the source space and the target space are both 2-dimensional, there is no need to train any autoencoder. A two dimensional extended OT map is directly computed. Our results are included in table 1, and benchmarks of previous methods are copied from Lin et al. (2018) and Xiao et al. (2018). Generally speaking, the samples generated by the proposed method can not only cover all of the modes, the quality of them is also better than others.
Besides, we experiment on stack MNIST dataset and CelebA dataset to further illustrate the performance of the proposed method, and the results are shown in Section C.1 and C.2 in the Appendix.

heading: 4.3 QUANTITATIVE COMPARISON WITH FID
text: FID is computed by: (1) extract the visual-meaningful features of both the generated and real images through the inception network, (2) fit the features in both the generated and real feature spaces with Gaussian distribution, and then (3) compute the distance between the two Gaussian distributions with the following formula FID = ‖µr − µg‖22 + Tr(Σr + Σg − 2(ΣrΣg) 1 2 ), where µr, µg mean the means of the real and generated features, Σr,Σg represent the variances of both distributions.
We report our results in Tab. (2), in which the compared data comes from Lucic et al. (2018)Hoshen & Malik (2019). In general, the proposed model achieves better than or comparable scores to other
state-of-the-art generative models. Theoretically, the FID scores of our proposed generative models should be close to that of the pre-trained autoencoders, and this is also validated in our experiments.
The autoencoder architecture we use here cannot find a good encoding for the CelebA dataset due to the limited capacity. But the FID score of the generation model is still approach to the autoencoder. In order to verify that with appropriate capacity of autoencoder, the proposed model works. We use the generator of DCGAN Radford et al. (2016) as the decoder of the autoencoder, then the reported FID score is 28.6, outperforming other models. Also, some of the generated images are displayed in Fig. 7 of the Appendix.
We also display the generating results for the four dataset in Fig. 5. It includes the original images, the best generating results of Lucic et al. (2018), including various GANs and VAE, the results of Hoshen & Malik (2019) and ours, row by row.

heading: 5 CONCLUSION
text: This work gives a theoretic explanation for mode collapse/mixture by Brenier’s theory and Figalli’s regularity theory of optimal transport maps. When the target measure has concave support, the OT map is discontinuous on the signular sets. But DNNs can only represent continuous functions, this conflict causes the both problems. In order to solve this problem, the AE-OT model is proposed by separating manifold embedding and measure transformation. The former step is computed using an autoencoder, the latter is carried out using the extended semi-discrete OT map based on GPUs. The model is tested thoroughly and extensively by both synthetic and real data sets. The experimental results validates the discontinuity of the OT maps and demonstrate the advantages comparing to the state-of-the-arts.

heading: Acknowledgements
text: The project is partially supported by the National Natural Science Foundation of China (61936002, 61772105, 61720106005), NSF CMMI-1762287 collaborative research: computational framework for designing conformal stretchable electronics, Ford URP topology optimization of cellular mesostructures’ nonlinear behaviors for crash safety, and NSF DMS-1737812 collaborative research: ATD: theory and algorithms for discrete curvatures on network data from human mobility and monitoring.

heading: A BRENIER’S THEORY
text: In this subsection, we briefly introduce the basic concepts and theorems in optimal transport theory, which comes from Brenier theory Villani (2008); Brenier (1987; 1991a), and discrete theory Gu et al. (2016).
Optimal transport Problem Suppose X,Y ⊂ Rd are two subsets of n-dimensional Euclidean space, µ, ν are two probability measured defined on X and Y respectively, with equal total measure, µ(X) = ν(Y ). A map T : X → Y is measure preserving, denoted as T#µ = ν, if for any measurable set B ⊂ Y , µ(T−1(B)) = ν(B). Given a cost function c(x, y) : X × Y → R≥0, indicating the cost of moving each unit mass from the source to the target, the total transport cost of the map T is defined to be ∫ X c(x, T (x))dµ(x).
The Monge’s problem of optimal transport arises from finding the measure-preserving map that minimizes the total transport cost.
(MP ) Wc(µ, ν) := min T#µ=ν ∫ X c(x, T (x))dµ(x). (2)
The solutions to the Monge’s problem is called the optimal transport map, whose total transport cost is called the Wasserstein distance between µ and ν, denoted asWc(µ, ν).
Brenier’s Approach Brenier Brenier (1987; 1991a) proved the following theorem: Theorem 1 (Brenier Brenier (1987; 1991b)). Suppose X and Y are the Euclidean space Rd and the transport cost is the quadratic Euclidean distance c(x, y) = 1/2‖x − y‖2. Furthermore µ is absolutely continuous and µ and ν have finite second order moments, then there exists a convex function u : X → R, the so-called Briener potential, its gradient map∇u gives the solution to the Monge’s problem. The Brenier potential is unique up to a constant.
Brenier’s polar factorization theorem claims that: for any measure preserving map T#µ = ν, T can be uniquely decomposes into the forms T = ∇u ◦ s, where s : X → X is a volume preserving map and ∇u is the optimal transport map under L2 cost. Therefore, the regularity of T can be determined by that of∇u.
Discrete Brenier’s Theorem Brenier theorem can be directly generalized to discrete target measure. Suppose the source measure µ is defined on a compact convex set Ω, the target measure ν =∑n i=1 νiδ(y − yi), µ(Ω) = ∑ i νi. The discrete Brenier potential is a piecewise linear function,
uh(x) = n
max i=1 {πh,i(x)} = n max i=1 {〈x, yi〉+ hi} . (3)
As shown in Fig. 3(a), the projection of the Brenier potential induces a cell decomposition of Ω, each cell Wi(h) := {p ∈ Ω|∇uh(p) = yi}, whose µ-measure is denoted as wi(h). Theorem 2 (Discrete Brenier Theorem Gu et al. (2016)). For any ν1, ν2, . . . , νn > 0 with ∑n i=1 νi = µ(Ω), there exists h = (h1, h2, . . . , hn) ∈ Rn, unique up to adding a constant (c, c, . . . , c), so that wi(h) = νi, for all i. The vector h is the unique minimum argument of the following convex energy
E(h) = ∫ h 0 n∑ i=1 wi(η)dηi − n∑ i=1 hiνi, (4)
defined on an open convex set H = {h ∈ Rn : wi(h) > 0, i = 1, 2, . . . , n}. Furthermore, ∇uh minimizes the quadratic cost ∫ Ω ‖x − T (x)‖2dµ(x) among all transport maps T#µ = ν. The gradient of above energy is given by ∇E(h) = (w1(h)− ν1, w2(h)− ν2, . . . , wn(h)− νn)T . The Hessian of the energy is given by
∂wi ∂hj = −µ(Wi ∩Wj) ‖yi − yj‖ , ∂wi ∂hi = ∑ j 6=i ∂wi ∂hj . (5)
The optimal transport map can be obtained by convex optimization. Furthermore, the optimization can be carried out using Newton’s method. The global linear convergence rate is guaranteed by the following theorem:
Theorem 3 (Kitagawa-Mérigot-Thibert Kitagawa et al. (2019)). Assume the cost function is quadratic distance, µ has convex support and also that (i) The probability density of µ is C0,α(Ω) for α in (0; 1]. (ii) µ has positive Poincaré-Wirtinger constant. Then the Newton algorithm for semi-discrete optimal transport converges globally with linear rate and locally with rate 1 + α.
Though in our method the gradient descend method is applied, the above theorem also ensures its convergence because of the convexity of the energy function E(h) we adopted.

heading: B FIGALLI’S THEORY
text: In this section, we show the fact that even for the case of single mode, the transport map may still be discontinuous, which will cause the instability of the training process of GANs. The arguments are mainly based on the regularity theory of transport maps developed by Figalli Chen & Figalli (2017); Figalli (2010) and so on.
According to Brenier’s Theorem Brenier (1987; 1991a), any transport map can be decomposed into a measure preserving map and a solution to the Monge-Ampére equation, which is the optimal transport map under the L2 cost function. Therefore, the continuity of the transport map can be reduced to the regularity (smoothness) of the solution to the Monge-Ampére equation. When the support of the target measure is convex and the density functions are smooth, Caffarelli showed the map is differentiable; otherwise if the target domain is not convex, Figalli showed the map is discontinuous, and gave precise description of the singularity set. In this section, we briefly introduce Figallis’ theory, and conduct an experiment using CelebA data set to show the existence of the singularity set, hence demonstrate the fact that the transport maps computed in GANs are discontinuous.
B.1 CONVEX DOMAINS - CAFFARELLI THEOREM
Let Ω and Λ are two bounded open sets in Rn, and let f : Rn → R and g : Rn → R be two positive functions such that f = 0 in Rn \ Ω, g = 0 in Rn \ Λ, and∫
Ω
f = ∫ Λ g = 1.
According to Brenier’s Theorem Brenier (1987; 1991a), there exists a globally Lipschitz convex function ϕ : Rn → R such that ∇ϕ#f = g and and ∇ϕ(x) ∈ Λ̄ for L2-a.e. x ∈ Rn. We say ϕ weakly solves the Monge-Ampére equation
det(D2ϕ) = f
g ◦ ∇ϕ in Rn, (6)
together with the boundary condition∇ϕ(Rn) ⊂ Λ̄. ϕ is called the Briener potential. As shown by Caffarelli [9], if Λ is convex, then ϕ is strictly convex, and it solves the Monge-Ampére equation 6. The regularity theory has been estabilished (see Caffarelli (1990a;b; 1991)), such as
1. if λ ≤ f, g ≤ 1/λ for some λ > 0, then ϕ ∈ C1,α(λ)loc (Ω).
2. if f ∈ Ck,αloc (Ω) and g ∈ C k,α loc (Λ), then ϕ ∈ C k+2,α loc (Ω), (k ≥ 0, α ∈ (0, 1)).
B.2 NON-CONVEX DOMAINS - FIGALLI THEOREM
However, if Λ is not convex, the regularity of the Brenier potential can not be guaranteed. Figalli gave examples in Figalli (2010), such that
1. Ω is convex, Λ is simply connected, but non-convex; 2. the density functions f and g are smooth, f ∈ C∞(Ω) and g ∈ C∞(Λ); 3. the Brenier potential ϕ 6∈ C1(Ω), the transport map∇ϕ is not continuous.
In this scenario, the transport map can not be learned using DNNs, and training process is unstable or the GAN model generates unrealistic samples.
Figalli’s construction Let ϕ : Rn → R be a convex function. Its subdifferential at a point x is defined by ∂ϕ(x) := {y ∈ Rn|ϕ(z) ≥ ϕ(x) + y · (z − x),∀z ∈ Rn}. ϕ is differentiable at a point x if and only if ∂ϕ(x) is a singleton. Figalli decomposes the set of non-differentiability points according to the dimension of the singular set:
Σk(ϕ) := {x ∈ Rn|dim(∂ϕ) = k}, k = 0, . . . , n. (7) For any k = 0, . . . , n, the set Σk(ϕ) is (n− k)-rectifiable. The set of reachable subgradients at x as
∇∗ϕ := {
lim k→+∞
∇ϕ(xk)|xk ∈ Σ0, xk → x } .
It is known that the convex hull of ∇∗ϕ(x), coincides with ∂ϕ(x). Theorem 4 (Figalli). Assume that there exists λ > 0 such that λ ≤ f ≤ 1/λ in Ω, λ ≤ g ≤ 1/λ in Λ, and that ∂Ω and ∂Λ are continuous. Then ϕ is strictly convex inside Ω. Moreover there exist two open sets Ω′ ⊂ Ω and Λ′ ⊂ Λ, with L2(Ω \ Ω′) = L2(Λ \ Λ′) = 0, such that ϕ ∈ C1,α(Ω′), ∇ϕ is a bi-Hölder homeomorphism between Ω′ and Λ′, and ϕ is an Alexandrov solution of 6 inside Ω′. In particular, Caffarelli’s regularity theory for strictly convex Alexandrov solutions of the Monge-Ampére equations applies to ϕ inside Ω′.
Figalli studies the singular set of ϕ in Ω, i.e. the set of points x ∈ Ω where ϕ is not differentiable, denoted as Sing. Figalli shows the singularity set has the following characterization,
Sing = {x ∈ Ω|∂ϕ(x) ∩ Λ = ∅,∇∗ϕ(x) ⊂ ∂Λ, ∂ϕ(x) 6⊂ Λ} . it can be decomposed into connected components Sing := ∪iSi. For planar case, Theorem 5 (Figalli Singularity Set). The number of connected components of Sing is at most countable. Moreover:
1. either Si coincides with an isolated point {xi} for some xi ∈ Ω, and in this case the boundary of ∂(xi) is entirely contained inside ∂Λ (so that ∂ϕ(xi) completely fills a hole in Λ);
2. or Si can be written as a disjoint union as follows: Si = ⋃ j γij ,
where γij : Iij → Sing are embedded Lipschitz curves parameterized by arc-lengh, Iij is an interval.
B.3 ELEMENTARY EXPERIMENTS
We have designed a numerical experiment to verify Figalli’s theorems in low dimensional case.
Fig. 6 shows another computational result, which demonstrates the singularity structure in Figalli’s theorem. The source domain Ω is the unit disk, the target domain Λ is with complicated geometry. The singularity set of the optimal transport map satisfies the description of Figalli’s theorem 4,
Σ1 = 3⋃ i=0 γk, Σ2 = 1⋃ j=0 xj .
∂ϕ(x0) fills the hole on Λ. For any interior point p ∈ γ1, ∂ϕ(p) is a line segment connecting two points on the boundary of Λ.

heading: C ADDITIONAL EXPERIMENTS
text: C.1 STACK MNIST EXPERIMENT
Experiments of varisou GAN models on stacked MNIST dataset are in consistent with Lin et al. (2018). For AE-OT model, we use the architecture shown in table 3 and table 4, with the decoder
architecture same as the consistent generator architectures in GANs, and encoder having mirrored architecture.
We test diversity of generated samples from our AE-OT method on stack MNIST dataset that consists of 128,000 samples in 1,000 modes with each sample stacking three handwritten digit images from MNIST dataset LeCun et al. (1998). Number of modes counts the amount of modes captured by samples produced a generative model. The reverse KL divergence is computed by first assign each samples to their nearest mode, and compute the KL divergence between histogram of sample count on each mode and the histogram of real data. We choose angle threshold θ̂ = 0.5 for AE-OT method. Details of network architectures are listed in supplementary materials. Experiments results are summarized in table 5, which show our method achieves best performance in terms of modes captured and reverse KL divergence on stacked MNIST dataset.
C.2 CELEBA EXPERIMENT
we evaluate our method on CelebA dataset by measuring collision probability in a batch of 1024 generated images of size 64-by-64. If a pair of identical images appear, a collision is declared, and thus higher collision probability means lower generation diversity. The same metric has been used in Lin et al. (2018) for evaluation of PacGAN. To make a consistent comparison, we design our autoencoder network with encoder having the same architecture as in previous work and decoder having a mirrored architecture of encoder. Angle threshold θ̂ is chosen to be 0.7 for AE-OT test. Results are listed in table 6, with corresponding images can be downloaded here. Results have shown that our method achieves best result in terms of probability of collision. Autoencoder network structures can be found at table 7 and 8.
C.3 LINEAR INTERPOLATION IN THE LATENT SPACE
Given any two images in the dataset, we can find the images between them by linear interpolation in the noise space because the one to one correspondence between µ mass centers in the noise space and the images in the dataset is provided by the proposed algorithm. For other generation models, though the interpolation can be done successfully in the noise space, they cannot find the correspondence from the noise space and the image space. The results of the linear interpolation are shown in Fig. 8.

heading: 25×273K 0.65 0.17 0

heading: 16×273K 0.86 0 0

heading: 4×273K 0.42 0 0

heading: 273K 1 0.33 0

Review summaries:
General Comments:  The generator in Generative Adversarial Networks (GANS) computes an optimal transportation from the noise distribution to the data distribution.  However, such maps are in general discontinuous.  Since deep neural networks can only represent continuous maps, this brings two problems: mode collapse and mode mixture. This paper approaches both problems using Figalli's regularity theory. They separate the manifold embedding (here an autoencoder maps input data to a latent space) from the optimal transportation (this map is found by convex optimization). Composing these two steps yields the proposed method. Their method basically avoids representing discontinuous maps by the generator. Empirically, the proposed method performs similar or better than state-of-the-art.